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FOREWORD 


This report presents the results of development and implementation of 
steepest-ascent optimization theory as applied to trajectory computation. 
This program depicts the results of the developments performed by 
efforts contracted by the Trajectory and Optimization Theory Branch of 
the Astrodynamics and Guidance Theory Division of the Aero -Astro dynamics 
Laboratory at MSFC. Questions and requests pertaining to this program 
should be addressed to Mr. Ron Toelle at the Trajectory and Optimization 
Theory Branch, Aero-Astrodynamics Laboratory, MSFC. 
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ABSTRACT 


LIFTING ROBOT is a minimum Hamiltonian- steepest ascent 
multistage lifting booster optimization program. It can simulate up to 
15 thrust or coast events, provides aerodynamic lift and drag as a func- 
tion of Mach number and angle of attack, and can maintain "g" limits 
by throttling. 

The payoff and terminal constraints can be selected from a 
library of 26 functions. In addition, intermediate point constraints, 
selected from the same library, may be imposed on the trajectory 
following any one of the thrust events. 

Through the use of input switches, a variety of vehicle parameters 
can be optimized in conjunction with the control variables \ -pitch and 
X - yaw. Tank limits of stages being optimized can be held and performance 
reserves as a function of AV, can be calculated. The impact point of any 
stage can be calculated and publishable tables can be printed. The working 
coordinate system and the environmental simulation conform to Apollo 
standards. 

This document contains the LIFTING ROBOT input description and 
an example problem. 
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1. INTRODUCTION 


The program described in this report is designed to optimize a 
large variety of multistage lifting booster trajectories. This objective 

* -T* 

is achieved through the use of the Min-H steepest-ascent trajectory 
optimization technique described in Reference 1. Briefly, the steepest- 
ascent technique requires that a reasonable, but nevertheless arbitrary 
choice of the controls be used to calculate a nominal trajectory. In 
general, neither the desired terminal state will result, nor will the per- 
formance index be optimum. Next, by solving the adjoint differential 
equations associated with the linearized perturbation equations about the 
nominal trajectory, impulse response functions may be determined for 
arbitrary small variations in the control variables, and influence coeffi- 
cients may be determined for arbitrary small variations in the control 
parameters. The choice of small changes in these controls, which 
simultaneously moves the terminal state closer to the desired terminal 
state and improves the performance index, 'is calculated. This change 
in the controls is added to the nominal control history and the process is 
repeated until the optimum is reached. 

The LIFTING ROBOT program can simulate a multistage lifting 
booster having up to 15 thrust events. The program can be used for both 
ground- launch and jump-start trajectories. Both the aerodynamics and 
throttling of the Space Shuttle are simulated. The working coordinate 
system and the environmental simulation conform to Apollo standards. 

Minimum Hamiltonian 
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The payoff and terminal constraints can be selected from a library 
of twenty- six functions. In addition intermediate point constraints, selected 
from the same library, may be imposed on the trajectory following any 
one of the thrust events. Also, the instantaneous impact point of the first 
stage may be forced to avoid specified regions. 

Through the use of input switches, a variety of vehicle parameters 
may be optimized in conjunction with the control variables pitch and 
X~yaw. Tank limits of stages being optimized can be held and performance 
reserves as a function .of AV, can be calculated. Staging on either fuel or 
time can be specified. Also, the impact point of any stage can be calcu- 
lated and publishable tables can be printed. 

The LIFTING ROBOT program has straightforward automatic 
convergence logic and a dynamic updating scheme for the control parameter 
weighting matrix. Convergence should be as reliable and sure as the 
original ROBOT program [9]. 

For the most part, this report is devoted to a description of the 
mathematical model used in formulating LIFTING ROBOT and to such a 
limited discussion of the logic structure as affords a complete description 
of program flexibility. 
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2. COORDINATE SYSTEMS 


The basic reference coordinate system in the LIFTING ROBOT 

AAA 

program is the inertial geocentric cartesian coordinate system X Y Z 

A 

shown in Fig. 1. This coordinate system has the Y axis pointing north, 

A A A 

the X and Z axes in the equatorial plane, and the Z axis in the meridian 
plane that contains the launch site at gyro release time. In the LIFTING 
ROBOT program gyro release time or guidance reference release (GRR) 
is a reference time occurring either prior to or at liftoff. 


Next described is the inertial cartesian plumbline coordinate 

AAA 

system x y z, in which the equations of motion are written. 


AAA' 

The plumbline coordinate system x y z , shown in Fig. 2 is 

AAA A 

formed from X Y Z by first rotating counterclockwise about X through 

• A 

0 1 and then clockwise about y through - 90. A z is the launch azi- 
muth angle and 6^ = tt/2 - where Q q is the geodetic latitude of the 
launch site. Both A^ and d Q are input quantities. 


The equations for transforming a vector from the X Y Z system 

AAA 

to the x y z system are 


— 


— — 

A 


A 

X 


X 

A 


A 

y 

= A 

Y 

A 


A 

Z 


z 
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where 


A = 


sin A 

z 

0 


cos.A 

z 


cos A sin 0 . 
z 1 

cos 0 1 

-sin A sin 0. 
z 1 


-cos A cos 0„ 
z 1 


sin 0^ 


sinA cos 0, 
z 1 


Since A^ is an input constant and 0^ is the complement of an input 
constant, the matrix A is also constant. 


In the plumbline system the position coordinates x, y, z and the 

AAA 

velocity components w, u, v are measured in the x, y, z directions, 
respectively. 


The plumbline system in LIFTING ROBOT differs from the 
Apollo 13 coordinate system [3] only in the names of the axes, i. e. , 


A 


~ A ' 

Z 


X 

A 


A 

X 


y 

A 


A 

. y_ 

Apollo 13 

Z 


L J ROBOT 


The third coordinate system used in LIFTING ROBOT is the geo- 

A A A 

centric spherical polar coordinate system 0 r 0 with coordinates 0, r, 

, AAA 

and 0. The 0 r 0 axes, shown in Fig. 3, point in the direction of increasing 
0, r, and 0, respectively, and may be formed by first rotating counter- 

A 

clockwise about Y through 0 and then rotating counterclockwise about 

A 

0 through 0. 
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AAA 

The projections of r on X Y Z are 


X 


sin 0 sin 0 

Y 

= r 

COS0 

Z 

_ 

sin 6 cos0 


and therefore the projections of r on x y z are 


X 


sin 0sin0 

y 

= r A 

cos 6 

z 


sin0cos0 


The transformation from x, y, z to ^ r, 6 is therefore 

-l/ a ll X+a 21 y+a 31 Z 


0 = tan 


a l3 X + a 23 y + a 33 Z 


------ 

= + y + z 


0 = cos ^(a^x + a y + a a)/r) 


where the a. . are elements of the A matrix described previously. 


AAA 

The equations for transforming a vector from the 0 r 9 system 

AAA 

to the X Y Z system are 


A 


A 

X 


0 

A 


A 

Y 

= B 

r 

A 


A 

z 


6 
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where 


COS 0 

sin0 sin 0 

sin0 cos 0 

0 

cos 6 

-sin 0 

-sin0 

cos0 sin0 

COS0 cos 0 


AAA 

and therefore the equations for transforming a vector from the 0 r 0 

A A A 

system to the xyz system may be written 


A 


A 

X 


0 

A 


A 

y 

= D 

r 

A 


A 

Z 


0 


.where D = A* B 



AAA 

Also, the inertial velocity components in the 0, r, 0 directions, 

u , v respectively, may be written 
s s 


w 

s 

11 

s 

V 

s 



If the multiplication A* B is performed and substitutions for 0 and 
0 are made in terms of x, y and z, the elements of the D matrix become 


d ll ’ (a 22 Z ' a 32 y,/r Sin0 
d 21 = <a 32 x " a l2 Z>/rsto9 

d 31 = <a 12 y-a 22 x)/r S me 

* T 

( ) Denotes matrix transpose. 
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12 


22 


32 


13 


23 


33 


x/r 

yfr 

z/r 

(d 12 cos 9 - a 12 )/sine 
(d 22 cos 6 “ a 22 )/sind 
(d 32 cos 9 ~ a 32 )/sine 
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3. GEOPHYSICAL PROPERTIES 



Described in this section are three geophysical properties of the 
earth which affect a rocket trajectory: gravitational acceleration, geo- 
metric form, atmospheric properties. 

3. 1 GRAVITATIONAL ACCELERATIONS 


The gravitational potential function, U(r, 0), used in LIFTING 
ROBOT is [4] 


ju 

U(r, e) = -f i + 


cj V 2 


+ 


3 

DJ 

35 


(l-3cos 2 0) +~ , 

5 \ r 


R \3 
e 


(3 "5 cos 0)cos0 


/K e\ 4 2 4 

— ■ (3 “30 cos 0 + 35 cos 0) 

\ r 


where CJ, H, DJ, R , n are input parameters which are, however, 

e e 

preset to 


CJ = 1. 62345 x 10 

H = -0. 575 x 10“ 5 
DJ = 0.7875 x 10~ 5 
R g = Earth equatorial radius 
= 6378165. m 

= Product of universal gravity constant 
and earth mass 
= 3. 986032 x lO 1 ^ m 2 /sec 2 
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The components of the gravitational acceleration vector in the 
plumbline system are calculated as the first partial derivatives of 
U(r, 6) with respect to the plumbline position coordinates, i. e. , g^, g^ 
and g are calculated as 




5U 


3r 

■JU 


30 

g x 


3x 


dx 



3x 

g y 


dU 

au 

Sr 


+ ^ 
ae 

ae 


^y 

3r 

s y 


a y 

g z 


SU 


or 



30 


3z 


3z 



az 


These equations may be rearranged into the form 


V 


X 


Pi ' 

t-* 

CO 

i 

g y 

= G u 

y 

- g to 

a 22 

_ g z_ 


z 


a 32 _ 


where 


G 

G 11 3 


1+CJ( 6 


/R \2 

V 


2 f^e\^ 2 

(l-5cos 0) +Hl— 7 I (3 -7 cos 6)cos9 


+ DJ 


R \4 


(— - (6 - 9 cos 2 0) cos 2 0) 


❖ 

These partials are given in Appendix B. 
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A* 


G, 


TO 2 
r 


/R \2 ,.K \ 3 „ 

2CJ cos0 - H (-“j ~ 3 cos 


R s 3 


f^e\^ 12 2 

+ Dj (— j (— ~ 4 cos 0) cos0 


* 0 ) 


Equations of the same general form are used in the Saturn V 
flight computer. [5] 


In the event that a spherical earth is to be simulated these equations 


become 


where 




X 

X 




’ G u 

y 

-V 


z 




and, of course. 


G TO 


0 


3. 2 GEOMETRIC FORM 

The earth is taken to be an ellipsoid, [6] as shown in Fig. 4. , 
which rotates about its polar axis with an angular velocity 
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Cbfl e 








The angular velocity, Q ,^ , and the flattening, f, are input constants that 
are preset to 


£2 = 7.2921158 x 10 ^ rad/ sec 

e 

f = 1/298.3 


The relationship between geocentric colatitude, 0, and geodetic 
latitude, 0 , is expressed by 

ctn 0 = (1 - f) 2 tan 0 

g 


The radius of the earth as a function of colatitude, R(0), is 

R(0) = (1 - f) R e / y , (l-f) 2 sin 2 0 + cos 2 0 


The derivative of R(0) with respect to 0, which is needed in order to 
calculate the time at which maximum dynamic pressure occurs, and 
altitude related terms in the adjoint equations is given by 


dR(0) 

d0 


3 

R(0) f(2-f)sin0 cos0 

(R ) 2 (l-f) 2 
0 


3. 3 ATMOSPHERIC PROPERTIES 


The earth is assumed to have an atmosphere which rotates with 
it at the same angular velocity, so that there is no wind over the earth’s 
rotating surface. 
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The PRA63 model atmosphere [7] routine on the MSFC system 
tape is presently used to calculate density, p, pressure, and speed 
of sound, s, as a function of the altitude, h, where h is calculated from 

h = r - R(0) 

During the adjoint integration, analytic derivatives of p, p , and s; 

and > respectively are calculated in PRB63 which differs 

p dh dh sdh 1 J 

from PRA63 only in that these derivatives are calculated. 


3-6 



4. CONTROL VARIABLES 


The time history of the orientation in space of the centerline, c, 

of the boost vehicle is determined by the control variable attitude angles 

X (chi-pitch) and x (chi-yaw) shown in Fig. 5. 

P y 


In addition to defining the position of the centerline, c, x and 

Pa a a 

X may be thought of as defining the auxiliary coordinate axes p c <fi 

y 

shown in Fig. 5. This auxiliary coordinate system is formed by rotating 

A A 

clockwise about z through Xp and then counterclockwise about p through 


A A A 

The equations for transforming a vector from the p c <p system 

..AAA 

into the x y z system are 


A 


A 

X 


P 

A 


A 

y 

= c 

C 

A 


A 

Z 


_ 


where 


C = 


cosy sinx cosx 
A P *P y 

-sin Xp cosXpCosx y 


sinx 


y 


■sinx p sinx y 
•cosx p sinx y 
cosx 
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The unit vector c which defines the centerline of the vehicle 
in the plumbline system is 


A 

c 


siny cosy 

A P A y 

cosy cosy 
A P 'y 


smy 


y 
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5. AERODYNAMIC FORCES 


The passage of the vehicle through the atmosphere gives rise to 
aerodynamic forces defined to act coincident with and normal to the 

A 

vehicle body axis c. In LIFTING ROBOT, the relative velocity, V , 

£ l 

is considered to be the velocity of the vehicle relative to the atmosphere. 


M 


l w \ 


jrQ,^sixi6\ 


w ' <a 22 Z ‘ a 32 y> °e 

u 

= 

u 

- D 

0 

= 

u ' (a 32 x ' a 12 z,Q e 

w 


w 


l ° ) 


7V (a 12 y ‘ a 22 X> \ 


It is convenient to define the magnitude of the relative velocity, V 


v H s I v H | =|A« ! +I 

A 

and a unit vector, V , in the direction of V as 

-K - x\ 





The aerodynamic normal force in LIFTING ROBOT is assumed 

A A 

to lie in the plane defined by c and V . This allows the aerodynamic 

-K 

forces to be described in the plumbline system using the nonorthogonal 

A A 

coordinates c and V and the total angle of attack, a. This situation is 
depicted in Fig. 6. 
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It is obvious from the figure that the vector of force due to F^^ 


is 


F. . = -F . . c 
AA AA 


and that the vector of force due to F^^. can be found by summing com- 

A A 

onents along c and V , i. e. , 

Jbi 

F AN = F AN (c °‘ “ q ' CSC “ V 


Also, a. is simply 


-1 A A 

a. = cos (c • V_,) 
K 


The dynamic pressure, q, is calculated as 


1 2 
q = 2 P V B 


The Mach number, M, is calculated as 


M = V R /s 


These can be used to calculate F. . and F ., T as 

AA AN 


F AA = qS C A (M,a) 


f an = qS c n ( m ' a) 
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where S is the reference area and C (M, a) and C (M, a) are assumed 

ii iN 

to be 

C = a(M) cosa + b(M) sin 2 a + c(M) 

A 

= a(M) sina - b(M) sina cosa 

This form for and is consistent with reasonable assumptions for 
lift and drag as is shown in Appendix D. The coefficients a, b and c can 
be fit directly to aerodynamic data for and or constructed from 
lift and drag data. 

Tables of a(M), b(M) and c(M) are provided in subroutine CACN 
in LIFTING ROBOT. These may be changed by the riser via a new data 
statement. 


As is shown below, the form of C . and CL,, is such that F „ 

AN AA 

and do not have to be constructed in order to write the aerodynamic 

force vector in the plumbline system. 


Summing F and F , the aero dynamic force vector, F , is 
AA AIN A 

found to be 



<f an cot “ ■ W ° 


- F AN csc a V R 
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Using the definition of and above gives 



F 


Ax 


Ay 


A v A 

= - qS(b + c) c + qS(b cosa - a)V 


R 


Az 
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6. BOOSTER CONFIGURATION 


In LIFTING ROBOT, the simulation of the thrust profile of a 
multistage booster is accomplished by synthesizing the profile from a 
sequence of up to 15 thrust events. By allowing the grouping joi these 
thrust events into stages to be specified byl input rather than by fixed 
internal logic, a great deal of generality is obtained. Fig. 7 depicts 
five thrust events as an example of such a sequence. 

The ith thrust event is characterized by five (seven in the 
atmosphere) items: 


1) F. 

i 

2) m. 


3)1 


‘'il 


V iZ 


v. 


i3 


V i4 


4) T. 

5) W J 

l 


Reference thrust per engine 
Flow rate per engine 

Number of inboard engines 

Cant angle of inboard engines 

Number of outboard engines 

W V 

Cant angle of outboard engines 
Thrust event duration 

Weight jettisoned at the end of each thrust event 


and in the atmosphere . . . . 


6 ) 

7) 



Engine exit area 
Aerodynamic reference area 






The convention used in LIFTING ROBOT for labeling thrust event 
and miscellaneous weight drop event times is also depicted in Fig, 7 . 

Thu st event times are labeled t^ and miscellaneous weight drop event 
times are labeled t„.. 

2i 

Note that there are six t . but only five thrust events of duration 
r The same is true of a picket fence, in that there is always one more 
picket than there are spaces. The T . may therefore be thought of as 
"spaces", and the i subscript of t^. as the "picket" number, with i = 1 
at the beginning of the first thrust event. 


From the figure it is apparent that the t^ are calculated as 




with t^ being defined as some input initial time. 

In addition to thrust event items. Fig. 7 also depicts two 
miscellaneous weight drops. A miscellaneous drop weight, as distin- 
guished from a jettison weight, can be dropped at any time. The ith 
miscellaneous weight drop is characterized by three items: 

1) W\ Miscellaneous weight dropped 

2) Time interval between beginning of nf th 
thrust event and miscellaneous weight drop. 

"W 

3) n. Weight drop time is calculated from the 

beginning of this thrust event. Can also be 
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thought of as "picket" number of the thrust 
event time to which is added to get 
miscellaneous weight drop time. 

The ith miscellaneous weight drop occurs at t^. The t^ are 
calculated as 


2i 


w 

= t, . + r . 


where 


j 


w 
n . 
1 


Note that with this definition, none, one or many miscellaneous weight 

drop events may be defined relative to any given thrust event, and may 

occur during that or any other thrust event. The only restriction being 

that t n . , , must be greater than t„.. 

2i + l 2i 


In LIFTING ROBOT the thrust events are grouped into stages 
through the use of the input array N0VENT. The first member of the 
N0VENT array should contain the number of thrust events in the atmos- 
phere, the second member contains the number in the second stage, etc. 
However they are grouped, all thrust events must be accounted fori 


6-4 



THRUST AND FLOW RATE 




6 . 1 


The use of the four numbers v^> to describe the 

effective number of engines leads to a rather cumbersome notation if 
they are used in each equation where the number of engines is required. 
Consequently, an effective number of engines operator, iu, is defined 
to be: 


v. 

1 


v. A cost* 

il i2 


+ 1 >. C cosv. . 
i3 i4 


V il 


+ V i3 


if v ^ multiplies or Ae^ 


• ■ i' 

if y multiplies nu or cm. 


The input thrust levels for first stage component rockets are 
considered to be nominal sea level thrusts. The total thrust, T, for 
all thrust events considered to be in the first stage is calculated from 


T = zu(F. + Ae.(p “ p )) s 
l l i s a 


T ~ v. Ae. p 
v i l a 


where p g is the sea level atmospheric pressure and 
thrust level. 


T 

v 


is the vacuum 


The input thrust levels for all thrust events other than those in 
the first stage are considered to be vacuum thrust levels, and the total 
thrust is calculated from 


cm. is defined in Section 7. 4. 2. 

i 
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T 


y.(F. - Ae.p ) = 
1 1 i a 


T - v. Ae. p 
v l l a 


while still in the atmosphere and 


T = n.F. = T 
11 v 

once the atmosphere is dropped. 

The total flow rate, m, in any thrust event is calculated from 


• * 

m = v. m. 

i l 


•6. 2 THRUST TABLES, FLOW-RATE TABLES AND 
DELTA-WEIGHT TABLES 

If the input variables JTHR^ = 0, thrust and flowrate of the ith 
thrust event are calculated as above. 

If JTHR. = +1. m. is calculated from a table of m. vs time and 
i i i 

F. is calculated from a table of F. vs. time. If JTHR. = -1, F. is 

i ill 

calculated from a table of F.. vs. time and m(t) is calculated as 
m(t) = m(t“At) +Am(t-At) - Am(t) 
where Am(t) is a table of mass, loss vs. time. 
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6.3 


THROTTLING 



If T /m in the ith thrust event is greater than GLIM = g • GLIMG. 

V 1 O 1 

where GLIMG is an input vector, the throttle value T is calculated as 

r = m • GLIM. /T 
l v 


and T and m become 
v 


T = T T = m GLIM, 
v v i 

m ~ t m = (m GLIM^/T^) m 
Throttling should only be used when JTHR. = 0. 

6. 4 ' STAGING ON FUEL 


If the input variable MSWCH. = +1, staging will be on time. If 

MSWCH. = “1, the ith thrust event will terminate when a weight loss equal 

to FUELG. has occurred. The actual burn time r. will then be calculated, 
l i 

FUELG is an input vector. 
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7. THE FORWARD TRAJECTORY 


This section contains the equations of motion integrated in 
LIFTING ROBOT, a description of the forward trajectory flight phases, 
the terminal functions which may be selected to define an optimization 
problem, and various users’ options associated with the forward trajectory. 


7. 1 THE EQUATIONS OF MOTION 


In general form, the equations of motion integrated in LIFTING 
ROBOT are written 


Vc 


Pa 


w = 


u = 


V = 


F /m + g 
x °x 

F /m + g 

y y 

F / m + g 
z to z 


w 


y = u 


z = v 




-m 


Impact point penalty function. Integrated if 

terminal constraint No. 26 is activated. See Section 7. 10. 


P n = T / m 
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T 1 A A 

*10 -m 11 -'" 7 ! 11 

P n = “(w-g x + u-g y + v.g z )/V I 

hi = ’¥f ■ 

p,„ = [(bcosa - a) (V -V„) - (b +c)(c -V )] 

13 m IK I 

Equations 9 through 13 are integrated only on a converged 
trajectory. 


The forcing functions F , F and F are 

& x y z 

F 

x 
F 

y 

F 

z_ 

with, of course, q set to zero when the atmosphere is dropped. 

Equations of motion 9, 10, 11 and 13 are easily derived by forming 
the dot product V -V and adding and subtracting T/m. p is known as 
the characteristic velocity, p^q is known as the turning loss, p^ is 
known as the gravity loss and p^ is known as the drag loss, p^ * s an 
aerodynamic heating indicator. 

The mass, m, is calculated from 


a „ a 

= [T - qS(b + c)]c + qS(bcosa - a)V 


R 


m = u + m 

a 
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where p is continuous and consists of the total propellants to be burned 
plus the payload, and 


J T 

m = (EW. + E W )/g 
a 1 . 3 fe o 

1 J 

The constant g Q relates mass to weight and is taken to be 


g Q = 9. SOSSSm/sec"' 


Since m is constant from one weight drop to the next, p = -m. for all 

tj!t u or V- 


By this artifice, the seventh state variable, p , is made to be 
continuous at all times, including those times at which mass is discon- 
tinuous. The primary advantage of integrating this particular choice 
of state variable is the ease with which mass can be reconstructed during 
the adjoint integration. Since p is continuous, it may be stored as a 
function of time on the forward trajectory, and therefore, even if m is 
a time varying function obtained from a thrust tape, the mass can be 
calculated during the adjoint integration by looking p up, updating m 
at the t and t^. and adding the two together. 

7. 2 GROUND-LAUNCH TRAJECTORY FLIGHT PHASES 


The flight profile of a ground-launch trajectory is separated into 
a number of phases. These phases are depicted graphically for a booster 
having n-1 thrust events. The symbols in Fig. 8 are discussed below. 
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7. 2. 1 GRR, At , t 
o o 


The input quantity At is the time interval between the time the 
coordinate systems are defined, GRR, and the lift-off time, t^. t is 
an input constant which is generally taken to be zero, t^, the time the 
first thrust .event begins, is set to t . If a non- zero value of At Q is used, 
the boost vehicle, which is fixed to the earth, will not be in the 
YZ plane at lift-off. 


7. 2. 2 Ground- Launch Initial Conditions 


The calculation of the initial, t , conditions for a ground-launch 

o 

trajectory proceeds directly from At and the input value of the geodetic 
latitude of the launch site, 6 . The geocentric colatitude of the launch 
site is 


0 = 7r/2 - tan *{(l-f)^tan0 ) 

o 

The radius of the launch site is R(0) and the initial velocity of the launch 
site is 


V = R(0)O sin 6 
o e 


The longitude angle subtended by the launch site ‘during the time 

interval At is 
o 


A0 = 0 At 
o e o 
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The initial plumbline velocity components are 


w 


cos A 0 

o 


o 

u 

= V A 

0 

0 

o 


V 


-sin A 0 

o 


o 


The initial plumbline position coordinates are 


X 


sin A 0 sin 0 

o 


o 


= R(0) A 

cos 6 

z 

o 


cosA 0 sin 0 
o ! 





The initial value of the seventh state variable is calculated from 

m , and the input value of initial mass, m , as 
a o 

p = m - m 
o o a 

The initial value of state variables 9 - 13 is of course zero. 

7. 2. 3 Lift-off -- Phase 1 


The interval t -» t , Phase 1 of Fig. 8, is the lift-off portion 
o Xj 

the trajectory. During this interval the control variables and x y 
chosen so that the launch vehicle will clear the launch tower. 


of 

are 
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Since the launch tower is constructed normal to the reference 
ellipsoid, the angular separation of the launch tower and north is 6 , 
where 


®l = e i 


The longitude of the launch site is 0 , where 

J-f 


0, = A0 +a (t-t ) 

L o e o 


A unit vector in the launch tower direction can be transformed 
into the plumbline system as 


- 


“ — 


- 

X 


sin0 sin0 

J-i J-J 


siny cosX 

*p y 

y 

= A 

cosd 


COSY cosX 

p y 

z 


cos 0 sinCL 
L L 


S in Xy 


/\ 

Therefore, in the interval t -* t T , c is calculated to be 

o L 

x 

a — 

c = y 

z 


Since the A matrix and 6 are constant and 8 depends only on 
a ,L ' ^ 

t; x > X an d hence c during lift-off are functions of time only, t is 
P 7 

an input constant. 
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7.2.4 Tilt-Over -- Phase 2 


During the interval t -» t , Phase 2 of Fig. 8, the vehicle is 

JLj _l 


A A 


caused to tilt over in the x y plane by calculating x and X . as 

p y 


X P = x (t 'V 

x y =° 


where X is a trajectory parameter and is an input constant. 


During X =0 flight, the equations given previously for c 


reduce to 


A 

c = 


sin X] 

cosX 

] 

0 


7.2.5 Pitch- Plane Gravity Turn -- Phase 3 

Following the tilt-over, a pitch-plane gravity turn can be flown 
in which 


X = 0 

y 

A A 

and Xp is chosen so that the angle of attack in the pitch (x y) plane is 
zero. This requires that during Phase 3 
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A 

C 


2 , 2 
w/y w + u 


u/|/ w 2 + u 


7.2.6 Optimal Flight -- Phase 4 

The pitch-plane gravity turn terminates at t , an input constant 

X 

marking the beginning of Phase 4. 

Prior to t the thrust vector control angles v and X are obtained 
X P y 

as a direct consequence of internal logic phases. After t , x^ and X 

are considered to be tabular functions of time. Time, X and X can be 

P y 

-specified at a maximum of 196 tabular points. These are broken up into 

four sets of control tables with a limit of 49 points each. Through input 

it is possible to specify the thrust event "picket" number at which control 

tables start and stop and the number of points in a table. Control tables 

should not continue across a coast or an intermediate point constraint . 

/ 

Since Simpson's rule is used to integrate products of impulse response 
functions during the adjoint solution, there should always be an odd number 
of points in a control table. 


The steepest ascent process converges on the optimal Xp> X 
time histories by updating the tabular control programs of Xp aa -d x^ 
(if specified by input) at each iteration. If the input quantity KWTA is 


Described in Section 7.3. 
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set to 3, both y and y are varied. If KWTA is input as 2, y is held 
*P ' A y y 

at zero and y is varied. 

P 


7.2.7 Exo- Atmospheric Flight 


At t the atmosphere is dropped and the 12th and 13th state 

variables are no longer integrated. The LIFTING ROBOT program shifts 

to a different set of derivative routines at this point in order to avoid 

bypassing terms that have to do with the atmosphere. The internal logic 

of LIFTING ROBOT is arranged so that t is a miscellaneous weight drop 

Q 

event time, i. e. , 



with I being the number of the miscellaneous weight drop even t which 
terminates Phase 4. The end of Phase 4 (or possibly Phase 3 - see below), 
marks the end of atmospheric flight and hence t^ must be defined on 
every ground launch trajectory. Note that this implies that there must 
always be at least one miscellaneous weight drop event. If none is 
actually desired, then a zero weight must be dropped. 


7. 2. 8 Elimination of Phase 3 or Phase 4 

If the input variable TCHFRZ is input equal to TTILT there will 
be no gravity turn and optimal flight will begin immediately after Phase 2. 
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If TCHFRZ is input greater than t„ there will be no optimal 
flight within the atmosphere and Phase 3 will last until the edge of the 
atmosphere at which point optimal control begins. 

7. 3 INTERMEDIATE AND/OR TERMINAL FUNCTIONS 

In order to define an optimization problem it is necessary to 
specify the trajectory constraints as well as the quantity to be maximized 
or minimized. Table 1 consists of a library of 26 (at present) inter- 
mediate and/or terminal functions and their formulas. Any one 
of these functions may be selected as the payoff and be maximized or 
minimized at the terminal time. Any physically realizable set of the 
remaining functions may be selected as trajectory constraints and imposed 
•at the terminal time.. In addition, any physically realizable set of these 
functions may be imposed as constraints at an intermediate time by 
inputting the number of the thrust event following which the constraints 
are Imposed as NVRST. Additional functions can be added by setting up 
a new GO TO number in ACSTOP and computing the magnitude of the 
function and its partial derivatives with respect to the first 7 plumbline 
states. 

7. 4 CONTROL PARAMETERS, PROPELLENT TANK LIMITS AND 

FLIGHT PERFORMANCE RESERVES 

In addition to optimizing the and X y time histories, the 
LIFTING ROBOT program can simultaneously optimize control- parameters 
selected by input from the control parameter library. 
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TABLE 1. FUNCTION LIBRARY 



Code 

No. 

Function Name 

Symbol 

F ormula 


Mass (Payload if 
payoff) 

m 


2 

Inertial Velocity 

V 

I 

,, 1 2 2 2 
V = V w + u -p v 

3 

Inertial Flight 
Path Angle 

t 


4 

Radius 

r 

r =; J X 2, ± y 2 + z 2 

5 

Energy 

s 

c - v 2 - 2p ~- 
3 I r 

6 

Angular 

Momentum 

C 1 

C — V A 2 4 - B 2 -4- C 2 where 

A ~ yv - uz 
B — zw - vx 
C ~ xu - wy 

7 

Inertial 

Longitude 

0 

. / a, ,x a y a z\ 

0 = tan ~ 1 ( 11 + 21 + 31 

Va^x+a^y + a^z/ 

8 

Inertial Heading 
Angle 

(3 

s 

9 

Colatitude 

e 

e=COa ' 1(a 12 X+a 22 y+a 32 Z)/r 

10 

Inclination 

i 

i = cos ^ (sin & sin # ) 
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Table 1 . { C ont 1 d) 



Code 

No. 

Function Name 

Symbol 

F ormula 

■ 

Line of Nodes 

CO 

OJ =. 0 -i-tan 1 (cos B tan 0 ) 

12 

Semi-Latus Rectum 

P 

2, 2 2 W 
p =r (w s + v s ) f 

13 

Eccentricity 

e 

e —J 1 -f- p. C 3 / ,U Q 

14 

Radius of Perigee 

r 

p 

These two constraints must be 

activated together. When they 

are specified by KCDPHI the 

state is propagated forward from 

tp until y = 0. The values of 

r at the points where y - 0 are 

compared and the smaller is set 

in r and the larger in T . A 
p a 

set of adjoint equations is inte- 
grated backwards from both 

r and r to t„. The values of 
p a f 

the X's at t^ are then stored into 

the vector partials of T and 

P 

r with respect to the state at 
a 

cutoff. A call to APPG initiates 
this process, and APPG is 
described in Appendix E. 

15 

Radius of Apogee 

r 

a 
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(3 


Table 1 (Cont'd) 


1 . 1 

Code 

No. Function Name 

Symbol 

Formula 

16 

True Anomaly' 

V 

q = (w g 2 + v g 2 ) * r / ju e 

_ i ( qu //w 2 -t-v 2 \ 
7] s= tan S s 3 J 

17 

Argument of Perigee 

a 

^ -if cos 6 \ 

a « Yj + tan ( 

\ sm y cos/3; 

18 

Outgoing Asymptote 

Asym 

A = F • s’ 

B — ,/r 2 - A 2 
C 1 

Asym = (Cj - B ^ 3 ) - r + A 

s is a unit vector in direction of 
asymptote 

19 

Asymptote Plane 

Asm PI 

Asm PI 3 C • s 

is the angular momentum 
vector 

20 

Rendezvous W 

Rend W 

Rend W = W-W . 

target 

21 

Rendezvous U 

Rend U 

Rend U = TJ-U. 

target 

22 

Rendezvous V 

Rend V 

Rend V = V-V^ 

target 

23 

Rendezvous X 

Rend X 

Rend X = X-X, 

target 

24 

Rendezvous Y 

Rend Y 

Rend Y = Y~Y, 

target 

25 

Rendezvous Z 

Rend Z 

Rend Z = Z~Z 

target 

26 

Impact Penalty 
Function 

PEN 

PEN = p 0 

O 
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7.4.1 Control Parameters 


n 

Jdr^| 


Table 2 contains the members of the control parameter library. 


Library No. 

Parameter Name 

Symbol 

1 

1st thrust event duration 

T 

11 

2 

2nd thrust event duration 

T 

12 

m 


* 

• 

* 

• 

* 

# 

• 

15 

15th thrust event duration 

T l, 15 

16 

Launch Weight 

m 

o 

17 

Tilt- over Chi' dot 

X 

18 

Launch Azimuth 

A 

• 


z 


TABLE 2 CONTROL PARAMETERS 


The library number of each parameter to be optimized is 
specified by putting a nonzero value into the equivalently numbered 
elements of the input array KDB. Thus it is the position of nonzero 
elements in KDB which indicates an active parameter. Although all 
T li ar e provided, a library number, only those r terminating after 
optimal control begins may be selected for optimization. 


7’. 4. 2 Propellant Tank Limits 


In a great number of real problems the total propellant in a 
given stage is fixed, albeit allocated among a number of different thrust 
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events- Also, since the available fuel and oxidizer will not, in general, 
be exhausted simultaneously when mixture.- ratio shifts are considered, 
tank limits in LIFTING ROBOT are based upon "critical" propellant 
rather than actual propellant. In what follows it is assumed that 
JTHR. = 0 for all thrust events involved, and that throttling does not 
occur. 


The r,. 
li 

relationship 


can be connected by logic so as to maintain the 


m = T. v . cm. r, . 
x .1 1 ll 

l 

where m . when tank limits alone are considered, is defined by the input 
x 

values of the critical flow rate cm. and the T , .. Since m cannot vary 

1 li x 

when the r^. are being varied by the steepest - ascent process 

S. v. cm. dr,. =0 
ii i li 

Therefore, all the connected thrust events cannot be optimized indepen- 
dently. One of the the jth, must be dependent and result from a 
choice of the others, i. e. , 

v. cm. 

dr.. = -£ — — dr 

^ i^j iLcm. 1 

J 1 3 

The procedure used in LIFTING ROBOT for specifying that the jth thrust 
event is connected to the ith with the ith being independent, i. e. , 
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KDB(i) f 0, is to put the difference between j and i into the same 
element of the input array KDT, i. e. , KDT(i) = j - i. One restriction 
on this procedure is that j must be greater than i. If no connection is 
desired, the appropriate element of KDT is set to zero. Note that if 

KDT(i) = j - i > 0, then KDB(j) must be zero since the same parameter 
cannot be specified as both dependent and independent. (If this requirement 
is not met, the program will print a warning, run a forward trajectory 
and go to the next case. ) If, for example, KDB(4) f 0 indicating that 
T ^ is to be optimized and if KDT(4) = 1, then is altered to keep 
m constant and KDB(5) must be zero. If, on the other hand, KDT{4) 

X 

= 2, then r is altered to keep m constant and KDB(6) must be zero. 

1 JU X 

If, however, KDT(4) = 0, then T ^ is optimized without regard to limits. 


As is implied by the equation for dT.^., the same thrust event 
can be specified as dependent by more than one independent parameter. 
For example, the input arrays 


KDT = 0, 0, 0, 4, 3, 2, 0, 0, 0, 0, 0, 0,0, 0, 0 
KDB = 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1 


indicate that m , v, r. t„ and t„ „ are to be optimized and that 
o A 14 15 16 


_ _ V m 4 V m 5 V m 6 

d 18 .. . d 14 . d 15 , d 16 


V m 8 


V m 8 


V m 8 
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It should be .noted that although the rationale for the development 
of the connection logic comes from the necessity of holding stage tank 
limits, the connection logic is independent of stage specification. 


7. 4. 3 Flight Performance Reserves 

Flight performance reserves (hereafter called FPR) is a name 
given to the propellants held in reserve on a design flight to provide an 
increment of velocity over and above the design velocity in the event it 
should be necessary on an actual flight. As such, FPR are jettisoned 
along with the jettison weight of the last thrust event and do not appear 
as part of the payload. Again it is assumed that JTHTh = 0 for all thrust 
events involved, and that throttling does not occur. 


The input quantity IPR is the number of the thrust event from 
which the FPR are withheld. If IPR = 0, FPR are not calculated. There 
are several accompanying requirements if IPR is not to be zero. First 
of all the IPR th thrust event must be in the last stage. Secondly, the 
maximum amount of critical propellant in the last stage, m , must 

X 

be input as WPMX. Thirdly, although the IPR th does not have to be 
the last thrust event, no thrust event which follows the IPR th may be 
optimized. 


The LIFTING ROBOT program calculates FPR on the basis of two 
input AY requirements. These are, AV^ to account for geometry per- 
turbations and AV to account for performance perturbations. FPR are 
P 

related to AV and AV through the equations 
g P 
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0 

-AV /Vex 

GPR = m (1 - e g ) 
c 

-AV /Vex 

PPR = (m - GPR)(1 - e P ) 
c 


FPR = GPR + PPR 


where m is the mass at cutoff of the IPR th thrust event, and 
c 

Vex = g I of the IPR th thrust event. Defining 
b o sp 


-AV /Vex 
k x = 1 - e g 

-AV /Vex 
* 2 ->-• p 


The FPR can be calculated as 


FPR = m k 
c 4 


where 


k 3 (1 “ k l* k 2 


Denoting IPR by j, and the mass at the beginning of the IPR th 
thrust event by m„ the cutoff mass, m^, can be written as 


m = m. - v . m. t, . 
c 1 1 3 1] 
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The problem of course is to find t ^ such that the sum of the critical 
propellant contained in the FPR and that consumed during the remainder 
of the last stage is equal to m . This may be written 

X 


m = E v. cm. t. + v . cm.(T .+ T ) 
. x . i i li 3 3 *3 P 


where the summation by i is over the thrust events in the last stage, and 
r is defined by 


r 

P 


FPR 

v.m. 
3 3 


m. 

= k 4<^- 
v m. 
3 3 



This leads to 


cm. 


’ = (m - k .m. J - E v. cm. tv .) 

li * / -i i \ x 4 i i i li 

J i;.cm.(l“kj J m. ifi J 

3 3 4 3 


If in addition to FPR, r ^ are optimized in the last stage, a 
different form of the 'equation for . is useful in the calculation of the 
steepest ascent influence coefficients. Denoting the mass at the beginning 
of the first thrust event in the last stage by m and noting that 


m. - m T 
3 L 


E 

i<3 


v . m. 
1 1 



d 

m 


where m is 
between 


the sum of the weights dropped (if any) in the interval 

and m ., the equation for r . may be written 
3 3 
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cm. ^ cm. 

T = (m - k ^-{m -m ) + £ (k. 1 m.~ cm.) y.r,. 

li . ,, . , x 4 . L 4 - i 

d y.cm. 1-kJ m. Ki . m. 

3 3 4 3 Ns 3 


i l li 


£ i/. cm. r,.) 
... i i li 
J> 3 


The situation' that exists when t,. in the last stage are optimized and 
FPR are calculated, and when there is KDT connection between the ith 
and IPR th thrust events is essentially the same, since is constant 
in either case. The difference is that in straight KDT connection the 
input thrust event durations define an m , whereas with FPR, m = WPMX 

X X 

defines r^. The similarity between IPR and KDT connection can readily 

be seen for the case where AV = AV = 0, in which case k . = 0, and 

gP 4 

for both IPR and KDT connection 

V . cm. 

dr . = - £ — L dr , . 

3 • • l 1 

J if 3 v.CToa. 

3 3 

The situations are in fact so similar logically that the LIFTING ROBOT 
program sets up and uses KDT connection logic whenever are optimized 
in a stage that has FPR. 


7. 5 JUMP START 

The input variable JUMP is the thrust event "picket" number at 
which a trajectory begins. If JUMP = 1, the trajectory will progress 
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through the ground-launch logic. If JUMP ^ 1, the trajectory will begin 
at that thrust event "picket" number. If JUMP < N0VENT (1) the trajec- 
tory will begin in the atmosphere. If not, the trajectory will begin out 
of the atmosphere. The starting state is specified through the input 
array VIV, the starting time by TZER0 and the starting weight by W01. 
When there is a jump start, t is set to TZER0 and all KDB and KDT 
below the jump start point are set to zero. IF VIY(7) = 0, the plumbline 
state w, u, v, x, y, z must read into VIV(l) - VIV(6). If VIV(7) = 2, 

V T , y, r 3 azimuth (A ), latitude (0') and oo must be read into 
I z 

VIV(l) - VIV(6). 

Setting 

a = 180 - A 

z 

and 

6 = 90 - 6' 

a) = tan ^(cos0tana) 

0 = 60 - ju 

Then, constructing a B matrix using 6 and 0 above and using the launch 
site A matrix, a D matrix can be constructed and used to calculate the 
initial plumbline state as 
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m 

|DRcj 


w 

o 


cosy cos a 

u 

o 

= V D V 

siny 

V 

o 


cosy sin a 

X 

o 


0 



= r D 

1 


z 

o 


0 



7. 6 10 KM, QMAX, 14 KM 


The program prints out as it crosses 10 km altitude, 14 km 
altitude and the point of maximum dynamic pressure. In order to find 
the latter, the time derivative of dynamic pressure, q, is used, q is 
calculated by forming the dot product of the partials of q wrt the plumb- 

line state, ^ , and the time derivatives of the plumbline state, p . 

SP 

That is. 


q 



where 
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pw 



PE 


oz 


p if (d 12'^f^ Ld 13 ) -P (a S2a- a 222> S2 e 


p a (d 22‘ d 23> ' p(a 12^‘ a 32— > “e 


H (d - dR . ( . e ) d ) - p( a w-a u) to 
p dh V 32 d6 33' 22— 12-' e 


and is calculated numerically using the PRA63 atmosphere routine. 
7.7 IMPACT POINT 


The LIFTING ROBOT program integrates the trajectory of the 

J 

jettison weight of the IMPth thrust event to impact (h = 0) if the 

input constant IMP is > 0. The forcing functions F^, F_^ and F^ on the 
impact trajectory are 


F 

x 

F 

y 

F 

z 


-pV R w 
- P V R n 

- p v rZ 
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where p = 0 for altitudes greater than 690 km and p calculated from 
PRA63 as a function of altitude for h< 690 km. Only the first 6 equations 
of motion are integrated on the impact trajectory. 

7. 8 ANALYTIC COAST 


Universal coast equations have been installed in LIFTING ROBOT. 
These equations solve the two-body problem in cartesian coordinates. 

A description of the method is contained in Appendix A. 


An analytic coast may be specified in any two zero-thrust thrust 
events. The method of specification is through NC0ST1 and NC0ST2. 

NC0ST1 = No. of thrust event for 1st coast 
NC($ST2 = No. of thrust event for 2nd coast 

If the analytic coast is used, the problem will print out START 
COAST and END COAST at the beginning and end of an analytic coast. 

No print-out will occur 'during an analytic coast. The analytic coast 
should be used with a spherical earth. It goes without saying that the 
analytic coast results in a great time savings if long coasts are considered. 

7. 9 RENDEZVOUS 

If terminal functions 20 through 25 are activated, a rendezvous 
will take place and additional input, which gives the position and velocity 
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of the target at t , is required. The position and velocity of the target 
in the reference coordinate system at the nodal point are calculated using 
the radius of the target orbit at the nodal point, its flight path angle, 
inclination, and velocity; RTGT, GAMTGT, INCTGT and VELTGT, 
respectively. 


z = RTGT 
N 

w N = VELTGT co s( GAMTGT) cos(INCTGT) 
u^ - VELTGT cos (GAMTGT) sin(INCTGT) 
v N = VELTGT sin( GAMTGT) 


The position and velocity of the target at any later time, t, are 
calculated by coasting t , seconds using the initial conditions specified 

laP 

above. 


The target is positioned in its orbit through the use of the input 
variable BTATGT. The input state defined above is caused to coast a 
length of time equal to 

BTATGT < RTGT 
7 tar t RAD * VELTGT 

where t is current time in LIFTING ROBOT and BTATGT is an angle 
turned- through in a circular orbit. 
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The final target state is then rotated into the plumbline system 
using the A matrix. The difference between this state and the vehicle 
state at t^ must be zero for a rendezvous to occur. 

7. 9. 1 Rendezvous Options 

If IAA is input as 1, the launch azimuth will be calculated internally 
to coincide with the azimuth of the target as it passes the launch site, i. e’. , 

A = sin 1 (cos (INCTGT)/sin 0. ) 
z a 

7. 10 THE IMPACT POINT PENALTY FUNCTION STATE VARIABLE 

If IPCNST f 0, the derivative for the eighth state variable is 
calculated by mapping the instantaneous state down to R^ using an 
analytic form for the f and g series. (See Appendix E). The elapsed 
time, DTI, to impact is also calculated. The latitude and longitude at 
impact ^ and 0^ respectively are then calculated. The calculations are 


- 


- 




X 

p 


X 



w 

y p 

= f 

y 

+ 

g 

u 

z 

p 


z 



V 


DTI = r S 1+% S 2+We S 3 

®p = sin' 1 ((x p a l2+ y p a 22 +z p a 32 )/R e ) 
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0 = tan 

P 


-1 


- x a 
P 


n +y - a 


x a. +y 


p 21 + Z p a 31 


- Reft+At + DTI) 


The following vectors are then formed 


A 0j 

« 

♦ 


0 - LL0, 

P _ 1 

• 

A0 W>I 


0 p ' LLC W 


and 


A 0, 


" 6 - LLA, 

1 


P 1 

• 

A 0 N0I 


6 - LLA AlrtT 
p N0I 


where LLA = TH1/RAD where THL is an input vector of latitudes and 
LL0 = PHI/RAD where PHI is an input vector of longitudes of vectors 
of impact area allipses. M?I is the number of active ellipses. 


There follows 


e. 

1 


= erf. A0 2 + p A0.A0. + cr 2 A0 2 , 

0i l 12i li 0i i 


i = 1, NQSI 
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2 2 

where or p 1ri . and are the arrays S2L$, RH12 and S2LA respectively 

01 12i 01 

which are calculated in AINIT. 


If e\ is greater than 1, e\ is redefined to be 


e . = 1. + e. (e. - 1) 

l k.' i 

-e . 1 

i 

z. = c e 
i n. 
i 


where c^ and are set to 1 and 2 respectively. Then 
i i 

N0I 

^ = £i h 

This has the effect of setting up a state variable which increases 

at two different rates depending on the closeness of the impact point to 

one of the impact area centers. The e are intended to be used to cause 

& • 
i 

very little penalty to be generated when far from, an impact area center, 

and the c are intended to weight the centers relative to each other, 
i 

The arrays S2L$, S2LA and RH12 are calculated as 


p^ = (LATWTH./LONWTH.) 2 - 1 

o: = R0TA/RAD 
(S. = RAD/LATWTH. 
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S2L0, = /? 2 (1 +p 2 cos 2 a.) 

ill i 

S2LA. = £ 2 {1 +p 2 sin 2 a.) 

l l l l 

RH12. = ^ 2 p 2 sin2a. 
i i i . i 

where LATWTH and LONWTH are input arrays of latitude width and 
longitude width assuming the impact area elipses were aligned with lines 
of constant latitude and longitude, and ROTA is an input array of ellipse 
rotation angles, 

7. 11 OUTPUT TABLES 

By inputting a nonzero value of NT ABLE output tables suitable 
for publication can be obtained. The output tables are printed only for 
converged trajectories. 

If tables are desired, additional input described in Appendix D 
is required. 
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8. THE BACKWARD TRAJECTORY 



Since the steepest ascent method converges on the optimum set of 
controls by adding beneficial changes to the nominal set, the effect of 
small changes in the controls on the terminal and intermediate functions 
must be calculated. This is accomplished through the use of the adjoint 
differential equations. One solution of the adjoint differential equations is 
required for each terminal or intermediate function being either optimized 
or constrained. The adjoint solutions proceed backward in time from the 
final time for the payoff and terminal constraints, and from the intermediate 
constraint time if there are intermediate constraints. 

The adjoint variables are used to form impulse response functions 

which give the effect of changes in x and X and influence coefficients 

P Y 

which give the effect of changes in the parameters. These impulse re- 
sponse functions and influence coefficients are then used in the steepest 
ascent formulae to calculate beneficial changes in the controls . 

Notation traditionally used to describe the adjoint solution is in- 
troduced below. 

0 The scalar payoff function. 

0 An m x 1 matrix of constraints. Includes both terminal 

and intermediate constraints. (Constraints satisfied when 

0 = 0 .) 

An m x 1 matrix of constant Langrange multipliers 

associated with the constraints. (This v should not be 

confused with the effective number of engines v defined 

i 

in Section 7 . ) 
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T 

<1? The augmented scalar payoff function 0+v 0. 

X A 7 x 1 matrix of particular adjoint solutions associated 

0 

with the payoff function. 

X^ A 7 x m matrix of particular adjoint solutions associated 
with the constraints. 

X A 7 x 1 matrix of adjoint solutions associated with the 

function <!?. When appearing without a subscript, X is 
the equivalent of the Euler -Lagrange variables used in the 
calculus of variations (c.o. v. X's) and are formed as 

X = X + X , v 

0 0 

8. 1 BOUNDARY CONDITIONS 

The boundary condition on the Euler -Lagrange variables X are 
known to be 


X 


T 


a$ 

ap 


Consequently, 


the boundary conditions on X^ 


and X, 


* 


are chosen to be 


a t = a0 i 
A 0 ap 1 1 = t 


and 



a£ 

ap 


t = 


[t for terminal constraints 

; 1 

'i 

^lj, j = NVRST + 1 for intermediate constraints 
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8. 2 THE ADJOINT DIFFERENTIAL EQUATIONS 


Defining the 8 x m + 1 matrix X to be X = [X , ! X, ], a vector 

z z T 0 • f 

variational Hamiltonian H can be defined as H = X p . The Euler- 

z z z x 

Lagrange or adjoint differential equations become 


SH 


X = -( — ") 
z dp 


Z\T 


(for backwards integration) 


m + 1 sets of adjoint equations are integrated backwards to t (one 

Q 

set for 0 and one set for each of them $’s). The reconstruction of the 

plumbline state, needed to calculate (dH /dp) during the adjoint run, is 

z 

accomplished by looking up stored values of the state as a function of time. 


The differential equations for the X are most conveniently written 

Zi 

in vector form. To that end, define 


w | 
X z j 



z 


* t~. in Section 8 is considered to be the time where optimal control 

begins. 
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where the superscript on X corresponds to the state derivative it multi - 
T 

plies in X p . The total X vector is now 
z z 


X = 
z 


m 


; x p 


Out of the atmosphere X^ can be written 


x v = x r 

z z 


x r = J* x v 

z z 


X m ={ 
z \ 


-[F F F ] • X V / m 2 
x y z z 


j-X m GLIM/T 
^ z v 


if not throttling 
if throttling 


XP= 0 

z 


See Appendix C for a definition of J. 
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In the atmosphere X^ 


can be written 



z 



“ a 2 V r 


+ a 0 X "t cl 
3 z 4 


A 

c 


* r , v , 3h X T X T , 

X = J X - a ( — ) “ 0 ! o (- 1 _ ) + a 

z z 1 dr 2 3r 3 


ct{ X 


v 


V 


3r 


+ a 


.a — .T 
3(c -V ) 
r 

3r 



f-[F F F ] • X V / m 2 if not throttling 
i x y z z 

[^2 ( ( (b+c) c - (b cos a-a) V r j* X^) - X^m GLIM/THH if throttling 



0 


where 


= H* + aa.((b + c)i^-M^± £ - ) i- ^)]<x v - 4) 

1 m dh m p dp 3M s dh z 


qS/r, >1 dp ,„,3b 3a ^ 1 ds 1n v A ^- 

- ^([b cos a - a) - ~ - M ( —..cos a - — ) - — J (A. V ) 
— p dh 3 M 3M s dh z R 


m 


O' 


2^ [(a ' M! !M COSa -iSa ))a i'^R ) + (2<b+C> + M4i |M 1)<jL I ‘ £)] 


a. 


& 

2 m 


(b cos a - a) 


0! , 


, psv, 

1 R , ^ V A , 

~ b (X • v„ ) 

2 m z R 


A A 

also, cos a = (c -V ), 

K 
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and 



and ^ 7 - are calculated by the cubic fit SPLINE routine, 
, dM SM 9M J 

J 777 and — ~ are calculated analytically by PRB 6 3 
dh q dh s dh 


In addition. 



{ £ 

oh) T = 

1 r 

< 

ar 

1 — 
i L 


l r 


if spherical earth 


a. 


* , 12 R(0) 3 

“ “ (r cos 6 - r a 2g .. ) — 

l 32 


f( 2 -f )cos 6 

R £ ')T" lf oblate earth 


av, 


(V, 


R 


) = £2 
R a |T e 


- a 22 - a 32 


w a, — v a 
— 32 “12 


- a i2 ~ - a 22 


z R 


a? 


= a 


.V . u 

X z a 22 X Z a 32 

. W .V 

X Z a 32 " X Z a i2 
X Z a 12 " X Z a 22 


a ( c •v r ) t 
ar 


; sin x a 
j y 22 


- cos x„ cos X,_ a 

P y oz 


■°e | sin x p cos X y a 32 ' sin X y a 12 

i cos x p cos x y a i2 ' Sin X p COS X y a 22 
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* 

Note that )P - 0 always, hence A* 5 is not integrated. If one ele - 

Zt z 

ment of z happens to represent the penalty function, then that element of 

A* 5 is 1 and all others are zero. However, if one element of z happens 

z v r ■ 

to represent the penalty function, the derivatives of A and A are af- 

p P ^ 

fected. This can be seen by noting that since A = 1, the variational Hamil- 
tonian will be 


H = A V • V T +A r " 
P P I P 


\ m • , . 

V X p m+P 8 


Therefore, to the standard derivatives of H with respect to the first seven 
states must be added 


3P Q * dp p rp 

[ <^r> : <t=> f 


These partial derivatives are calculated analytically in IMPPT. (See 
Appendix F). 


8. 3 IMPULSE RESPONSE FUNCTIONS AND I INTEGRALS 


The impulse response functions for x and X axe defined by 

p y 

the equations 


T 

G = — 

zp dX p 

T hiX z ? ) 

G = 5 — 

zy dX 


T (r W cosx cos x -r^cosx sinx ) 

m z A y p z y P 


T (-r W sinx sinx - T U sinx cos x +r V cosx ) 

m z A y p z y p z y 
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where 


T 
= ( m 


m 


out of the atmosphere 
T -qS (b+c) 


m 


in the atmosphere 



xT out of the atmosphere 

A V + (\ V ' V_ ) V— in the atmosphere 

z z R T R 
m 


Impulse response functions are calculated at every tabular point in use in 

the X “X control tables. 

P y 


Denoting G by 
z 


G = 


G 


0P 


G 


0P 


if KWTA = 2 


G‘ 


rv 

0P 

G 

0y 

G 

! 0 P 

L. 

G, 

0y 


if KWTA = 3 


The m + 1 x m + 1 matrix of control variable "i" integrals is 
calculated during the backward trajectory as 


•' j a * I a > 

00 ; 00 ! tp rp i a 

= f 1 G W G dt, I (tJ = 0 
J z a z zz ± 


zz a , a 

. 00 ! 00 


Q 


where W 1 is a time varying weighting matrix defined in Section 8. 6. 

cl 
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If there are intermediate constraints, the above definitions of 

cl 

I may be used provided that after the intermediate constraint time the 
z z 

elements of G corresponding to the intermediate constraints are taken 
z 

to be zero. 


8. 4 INFLUENCE COEFFICIENTS 


The influence coefficients for a parameter give the changes in 
trajectory functions resulting from a unit change in that parameter and 
hence may be considered trajectory to trajectory partial derivatives. 

The influence coefficients for lift-off weight and tilt-over x are calculated 
using numerical derivatives; whereas, those for the are calculated 

using analytic partials. 


8.4,1 Influence Coefficients for Lift-Off Weight, Tilt-Over x and A 

z 


If the launch weight is to be optimized,, 

from t -* t with m changed by ±im.. 
o Q o u 

for the 16th parameter are then calculated as 


two trajectories are run 
The influence coefficients 


L 

zm 


p + (t o> - p'V 


2 A 


mo 


X (L) 
z Q 


where p and 
tive variations 


p refer to the plumb line state from positive and nega- 
of Am^ respectively. 
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If the tilt-over x is to be optimized, two trajectories are run 

from t t with' X changed by ± Ax . The influence coefficients for 
o y 

the 17th parameter are then calculated as 



+ 


p V ~ £ 

2 Ay 


X <t_) 
z Q 


If the launch azimuth A is to be optimized, one trajectory is 

z 

run from t to t„ with A changed by AA . The influence coefficients 
of z z 

for the 18th parameter are then calculated as 


zA 


■[(-s-r-m 


/A A 


8. 4. 2 Influence Coefficients for the r, . 

li 


The calculation of the influence coefficients for the t^. proceeds 
through three phases. In the first phase the influence coefficients for the 
effect of shifting the time at which a discontinuity occurs are calculated 
as 


L 

21 


T 

A p A 

z 


for the effect of shifting each and as 


. T 

Y . = Ap A 
Z] z 
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H 

iDftCj 


for the effect of shifting each t„. where Ap = [p - p + ] is the discon- 
tinuity in the plumbline state derivatives resulting from a discontinuous 
change in either thrust or mass or both. The Ap are calculated and 
stored during the forward trajectory; and L^. and Y^_. are calculated 
and stored during the backward trajectory. 


When, t . is the final time, p is set to zero. When t, . is the 
lx + ii 

intermediate orbit time, p is set to zero for the multiplication of those 

columns of \ corresponding to the intermediate constraints. In the 
z 

second phase cognizance is taken of the fact that the t are pinned to 

W n W J 

the t„ via r . . Since the r. are constant, 
li ] 3 


dt = dt H , i = NOWD(j) 


and therefore the following additions are performed in sequence with j 
running from l-*nw 


L . = L . + Y i = NOWD(j) 
zi zx zj 

where nw is the total number of miscellaneous weight drop events. 

In the third phase cognizance is taken of the fact that for the 
to be parameters 



Sr. dr. 

i 1 


1, 


3 > i 


and therefore the following additions are performed in sequence with i 
running from nv-1 1 
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L. . = L/ . + L . , j = i + 1 

Zl Zl ZJ 


where nv is the total number of thrust events. 


8. 4. 3 Influence Coefficients with Tank Limits and FPR 


If flight performance reserves are withheld from the jth thrust 
event, the influence coefficients for launch weight and r^. become 


L = L - v.rh. ~ — - — r L . 
zm zm j j (l-k 4 ) z 3 


v. cm.k 

L ’ - L . + 1 --- L . (i < L ) 

zi zi v.m.(l-k.) zi L 

3D 4 


where i is the first thrust event in the last stage, and 
JLj 


cm.m. v.cm. 
L . = L . - (1 - k . ■ - 1 . 1 ) r 1 


zi zi 


4 crn.m v.cinil - k ) zj L 
i 3 i J 4 


L ,(i r * i< 3> 


In the above equations, when and if the element corresponding to 
payload is augmented, k^ is set to zero. 


The proper augmentation of L . ,when tank limits alone are con- 
sidered and the jth thrust event is connected to the ith via KDT, is 


v .cin. 

L . = L . - L . (i< j) 

zi zi v.cm. zi 
3 3 
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Note that this result can be obtained from the last equation given above for 

FPR if k, and i r are taken to be zero. 

4 L 

8. 4. 4 Influence Coefficients with Staging on Fuel 

If any of the thrust events terminate on fuel. i. e. , MSWCH < 0, 

then the launch weight influence coefficient and the L . must be altered. 

zi 

This is becuase a change in mass causes a change in burn times (if there 
is throttling). Since fuel is constant, the difference between initial mass 
m. and final mass m^ of the thrust event is constant, i. e. , 

fuel = m. - m„ = const 
i f 


Therefore 


dm. = dm, 
i f 

Also, the mass when throttling begins is always constant since throttling 
begins when T/m = GLIM. Thus 

m. - m ( t - t.) 
l T l 


dm. 

l 

m 


and hence 


T/GLIM = 


dt 
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Since dm/dt when there is throttling is 


dm . GLIM 

— - — = - m — — — m 

dt T 

v 


v , m ( GLIM) „ , 

m, = Ttrrrr exp(- (t e - t_)) 


f GLIM 


v 


“f T 


and 


Therefore 


, , m GLIM _ , „ dm i v 

dm f = dm. m. f (dt ) 

1 v m 


dm. 

1 


1_ 

m 


m GLIM ni 


m. 

i 




Notice that if there is no throttling 


dt. 


- — = 0 since m. = m. 
dm. i f 

l 


If the launch weight parameter is active and there is staging on 
fuel in say the jth thrust events 


L = L + E L . 
zm zm . z~i 
J J 


dt. 


\ d 4i 


a 
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In addition, the following is done for each active tv incluence coefficient 


L ■ = L + E L 
z. z. . z. 

1 1 J>1 3 


f dt i 

(m ). ( - — 
z. f i Vdm. 


dm./j 


8. 4- 5 Parameter I Matrices 


Grouping the active parameters into a matrix = [L^ < L^] 

the m + lxm + 1 parameter "i" matrix can be formed as 

zz 


zz 


T b 

_b 

*00 

1 T 

00 

_b 

■ T b 

•©- 

1 

W 


T -1 
= L W, Lr 
z b z 


•where is a weighting matrix defined in Section 8. 6 


8. 5 STEEPEST ASCENT FORMULAE 


Denoting the vector of active control parameters by b, and the 
vector of active control variables by a, (if KWTA -2, a is a scalar 
equal to X p ) the steepest ascent formulae for the changes in the controls 
are 


6a 


db 


+ W * (G, - 
— a 0 




^0 ^00 ^00 
^“0 ^00 *00 


> E - w a ' 

> E - V 


c-r 1 

V# 

L 0 1 00 1 


k0 

k0 
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where 


*00 
X 00 

In the control equations above the plus sign is used when 0 is to be 
maximized, the minus sign is used when 0 is to be minimized, 0 < E < 1 
is a constant chosen to aid convergence, 0 is the vector of terminal 
constraints violations, and k is the decimal fraction of the constraint 
violation to remove. 


i a + i 

00 # 

i a + 

00 00 


If there are connected thrust events involving tank limits only, 
the dr^ for optimized thrust events will appear as elements of the db 
vector and the corresponding dT^. must be calculated as indicated in 
Section 7. 4. 2. 


The changes in the controls calculated using the above equations 
are then added to the nominal set to get the controls for the next iteration. 


The change in the payoff function 0 resulting from the control 
changes is 


“ ± ^ 00 “ *00 ^00 \< j ) E I 00 X 00 
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where the sign is chosen as before and 





8. 6 THE AUTOMATIC CONVERGENCE SCHEME 

It is the function of the automatic scheme to pick k, E, W ^ and 

-1 a 

in order to speed convergence* and to terminate a run when it does 

converge. The logic for picking k and E is straightforward and 

is directly related to iteration number. On the first iteration, E is 

set to zero and k is chosen such that .25< k< 1. A starting value of 

k can be input as DP2, however the program will ignore k<.25 or k> 1, 

3{C 

If k is input > 1 the iteration number is advanced to 2. On the second 
iteration k = 1 and E = 0. On the third iteration k = 1 and E = QY/2, 
where QY is an input constant which should be chosen 0 < QY < 1. On 
the fourth and subsequent iterations k = 1 and E = QY. 


The choice of the weighting matrices wj' and is also 
dependent on iteration number. On the first iteration W^ 1 is chosen to 
be 





m 


T 


m 

T 

0 

0 

m 

T 


if KWTA = 2 


if KWTA = 3 


If k starts at . 25 the k sequence is . 25, . 5, 1. ... with E = 0. , 0. , 0. , 

QY/2, QY . . . 

sL- Ox 

m/T is defined to be 1 / T 

m 
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and W, is chosen to be 
b 


W. 


-1 

b 


W 1 P 1 0 0 


W 2 P 2 


o w p 

n n 


where the W., are an input set of weighting numbers for the np active 
parameters. The W are input as WIBT and should generally be left 
at their preset value of 1 unless experience dictates otherwise. On the 
first and subsequent iterations the P. are chosen automatically so that 

* cl 

the largest contribution of the ith parameters to the diagonal of is 

equal to one. Denoting the influence coefficients of the ith parameter 

on the constraints by L^ , the P.th scale factor is 

J ij) i 


Max 

j = 1, m| 


Vi 


For the second and third iterations the constant Lagrange multipliers 
on the constraints, V, are 


= “1 
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thereafter, v is formed as 


v " \<p + * 

where the minus sign is used of maximizing and the plus sign if minimizing. 

Once v has been calculated, min-H on the control variables can 
begin since the Euler-Lagrange multipliers X can be formed as 

x - VV 

the variational Hamiltonian H has 

T 

H = X p 


for the first partial of H with respect to Xp and X as 


5H BH 

% ax 


and the second partial of H with respect to 



and x 

y 


as 


H 

aa 


S 2 H 


9X 


P 


5 2 H 

s V* y 


3 2 H 


9x p 5x y 


5 2 H 

^x„ 2 
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Therefore on the second and subsequent iterations is taken to be 

W" 1 =+ H~ X 


a aa 

where the minus sign is used of maximizing and the plus sign is used 
minimizing.. 

The elements of H are 


if 


9 2 H 


ax 


p 


dH 


ax ax 

p y 


a 2 H 


ax 


aa 

_w 11 

= -T (T cosx smx + r cosx cosx ) 

m y p A y A p 

rr-iW __H 

= -T (T sinx cosx “ F sinx sinx ) 
m A y A p y P 

= -T (T7cosX sinx H-r^cosX cosx +r V sinx ) 
m 1 A y A p y A p A y 


where 




pW 


if 



0 


0 

r U 


pU 

+ 

if 



0 


0 

r v 


r v 


if 

- **- 


0 


_ 0 _ 
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If H is ill conditioned, X and X satisfying H = 0 are used 
aa P y a 

to calculate a backup H having elements 

ctcl 


ax 2 


= + T m [(T w ) 2 +(r u ) 2 ]/y(r w ) 2 +(r u ) 2 +(r v ) 2 


d 2 H 

S V X y 


= o 


= + T ■\j^f+(T a ) 2 +(r v y 

2 m v 




where the minus sign is used if maximizing and the plus sign is used if 
minimizing. 


If KWTA = 2, H is 
aa 


H = -T <T W sinX +r Ll cosx > 

aa m p p 


-U 


with backup 


h = +t Vcr w ) 2 + <rV 

aa m V 

If a backup H matrix is used, the output quantity KAT will be 
1; otherwise KAT = 0. KAT should never be 1 on a converged run. 
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On each iteration a normalized total influence coefficient for 
each parameter, L 1 is formed as 




Prior to the fourth iteration the input values of W. are used in 

-1 1 
the construction of W. . For the fourth and subsequent iterations each 

b 

W. is altered according to the following logic: 


unchanged if 
W\ unchanged if 


L 1 ] < . 005 


present 


last 


otherwise 


W. = 2W. if [L 1 - L X t ^.1 <i L 1 ! 

i i ~~ present — Last — 


W i ■ w i /2 if ^ — ^present” —'last 1 — ^ ^ 


present 


present 


The are printed out as WIBT between iterations. 

This dynamic updating of will generally insure smooth 
convergence of the parameters. The relative magnitude of the W\ on a 
converged run can be used as a guide in picking input WIBT. 
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If there have been at least 3 iterations, if | dm Q I <100 kg, if 
I | < . 00002 radians, if all j dr^. J < . 5 seconds and if in addition all 
J L 1 ! < . 005, the parameters are considered to be converged and the output 
quantity BETCON will be T; otherwise BETCON will be F. 


The convergence test for the control variables x a nd x Is 



max 


d X, 


max 


Max 
overall 
points in 
chi-tables 


H _1 H j < . 005 
aa a 1 


This implies that the max deviation of either x^ or x from the 
optimum anywhere along the trajectory is less than . 005 radians. The 
max deviation in X from the optimum is labeled DEL CHIP MAX in the 
output and the max deviation in x is labeled DEL CHIY MAX. 

y 

As soon as I dx ! < . 005, I dx | < * 005 and BETCON 

1 p 1 max 1 v 1 max 

is T, a run is considered converged. A final forward trajectory is 
then run at the input print interval, integrated impact (if any) and output 

tables (if any) are run from this trajectory and then LIFTING ROBOT 
looks for input for the next case. 
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APPENDIX A 



UNIVERSAL COAST EQUATIONS AND THEIR 
APPLICATION TO TRAJECTORY OPTIMIZATION 


First, we begin the analysis of the universal coast equations 
with a statement of plausibility: Since any conic is completely described 
in terms of six elements, generally c , c , etc. , it is not^ unreasonable 

JL O 

to search for a solution of the equations of motion which has as the six 

elements the state vector at the initiation of coast, i. e. , x , x . 

o o 

A. 1 EQUATIONS OF MOTION (SCALAR FORM) 


The equations of motion in terms of the radius vector x are: 


x 


JUX 

3 

r 


( 1 ) 


where 


r = (x • x) 


1/2 


Calculating the second derivative of r from Eq. (2) gives 


rr = x * x 


and 


.. If. 2 2 u, , 2 , ul 
r = — x — ^ - r + ^ 
r L r r -J 


r >- r 

Since the energy a = , Eq. (4) becomes 


• • 1 / .2 , H x 

r = — (o' - r + ^ ). 

r r 


( 2 ) 

(3) 

(4) 

(5) 
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Now make a change in the independent variable defined by the equation 


dt - r d "'f 


(• 6 ) 


Consequently, 


dr _ d_r drt _ Y 

dt ~ dt av- " 

d Z r _ d (r r ) d_t „ Vy-Z+y- 2 *- 

TF ' dt d> 

If we denote differentiation by ^ with a prime, then, from Eq. (7) 


( 7 ) 

( 8 ) 


r ~ r / 


(9) 


and hence 


r'‘ - y r 2 t ^ 


(10) 


By substituting Eqs. (5) and (9) into Eq. (10), 


y" =. r<A - x.. - + /i t 1 
r r 


/Z 


(H) 


or 


r = r<* i- M 


(32) 


Now we seek an expansion of r in terms of 7 * about the point 


as 


^ +• r 0 ' v- * 


r o" y z * v" y 3 + r 0 ,Ji ‘ v 4 


■H - 


(13) 


Zl 


31 


We already have Y 0 " and the higher derivatives of r c are given by 
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Therefore Y" 

r = r 0 +- 
+ 


r'" - <*r' , 

y = cAr" - <x z r r <A JJ. , 

y‘"“- eA z r'j y 

r"""= <* z r" = <X 3 r 4 <* 2 >u , 

etc. 

may be written 


V? + 


V ^ 


4-/ 



^ 

n 

zl 

Zl 

3> 



» 4l 

f r 0 S t +• 

>■€}<} 


The terms in Eq. (15) may be grouped in the following way: 




t* V a 




r -- r Q 

i -t- 

1 4 - 

2 ! 

4 1 * 

* » * • <* « 

*>•' J 



4 t 

ca4 3 

**+*' 

i 



3! + 

S! + 

7 / 4 . * - 


t M 

' if, 

L zf 

*4* h 

4! 

+ 

^ ! 

» « •* 


Define : 





> 

5 0 = 

i c 

* V" , 
2 ! + 

i-* ^ 
— +r- * 

; - ■— / 000 

6 l 


7 / 

* ■/-* 


^ f 7 


$ ~ 

r f 

3! + 

St + 

yr 4- '« * 

> 

= 

m 

t «‘l f * 

+ 

-p * C c* 


2 ! 

4-1 

(ol 



■Sjj - 

p 

f- 



31 

Sf 

7f 




(14) 


(15) 


(16) 


(17) 
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In fact, Eq. (17) can be written in general form as 



y^ pt ~ n 

(2> p r n 'i l } 


( 18 ) 


where the number of terms desired, depends on the number of 

significant figures desired. 

Using the definitions in Eq. (17) and defining (T Q = T 0 - T 0 T 0 , 
Eq. (16) may be written: 


Also, since 


r s r Q S c t (T 0 S i i- jjlS z (!9) 

i --J r<jp , 

t = -t 0 t r 0 S, i- 62 +■ M<S - 3 - ( 2 °) 

Eq. (20) is solved iteratively for that value of J f which yields the desired 

d 4' 

coast time (i - t 0 ). The fact that -y j s useful here. Only the 
highest two values of need be calculated as a series since 


S n ~ $ n -t 


Z * 


n * 


(2D 


A. 2 EQUATIONS OF MOTION (VECTOR FORM) 


Now go back to Eq. (2) and proceed to derive an expansion for 
the state vector X. , X similar to that derived for Y . Again we note that 


and 


dx dx di c t 

d 7 " dT dY" Xf 


( 22 ) 


d\ _ 
d V 2 


d 

d-t 


Cxr) 


57 


“2 * » >< 
x ( t X r r = X - 


(23) 
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Ther efor e. 


x'r' 


X // _ - ^ * +- 

a .. r r r 


(24) 


r (- + x y/ r ' + r"x ') - (x'r ' - mx)y' 


/// _ - xayk' t rr' -prx / (rtx-j-/A) - fX / r / - / Mx)r' / 


1 1 > 


& , P, , & , , % P .P, / 

~ /xrx +- x / r - /U\r + xjrjxrx/JL f - Xrr t/AKf 


€ / 


X 
X ' 
X' 

etc. 


.- nt 


X <X 


/ . t 


X 1 (A = <* <x , 


r 


a' m , 


(25) 


(26) 


Expand X. and X in a Taylor series about the initial point: 


X - X 0 t- X, 


/ , x " tf 2 - x 5^ 3 

, f t *° X-f %> ■ X - -t- * . ' 

> / 2l 3 / 


r / // / X r ^ ^ 

X --K 0 tX^(- - 4 T7 2 - f 


//✓✓ 


P, 


l (27) 


2 L_f , 

£/ 


J 


Substituting Eqs. (25) and (26) into X we obtain 

A - *« + X.'f * *&- tr- ^ 


r„ 


2 ' 


f x 


/ <A f 


3 


XoV/ A _ ^X^ A P 

r c 41 


, w . ^ - » , A — ' *C ~ ' J_ 

37“ * 1T“ "4T ^ 


(28) 


Collecting terms, Eq. (28) becomes 

M fV . 


X 


*•['- P ! zr^^tr + • 


•)] 




( 29 ) 
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or 


* - *o ( 1 ~ ) f r c t rj $ z )' (30) 

r 0 f o 

Since 

&~0 ' > 

X 0 j 

^ i i~ c — £ ■* Y. 0 " yUvSj ., 



the equations for the state vector x may be written 

X- fx o .+ JK » 

where 

■f - I - M ^Z 
To 

and 


(31) 


J 


3 r k~i 0 ~ M-S 3 * 

Substitution of the expressions for the higher derivatives of X into the 
expansion for X results in ; 

mSi 


(32) 


x %rx = (- f (r- M 5 z )K 0 


Hence 


where 


r 0 

X - -Px 


o fcjx, 


(33) 


(34) 


b df M^i 

dY d-fc rr c 

o = ^ at , j _ mS 2 
^ d Y 3T 


"T 


(35) 


We have succeeded in the first of our objectives, namely, to express 
the state at one time as an explicit function of the state at another time. 


A- 6 



A. 3 TRANSFORMATION OF ADJOINTS 

We proceed now with, the second of our objectives: to develop 
the partials of the state at one time with respect to the state at another 
time. This is necessary to transform the adjoints across the coast 
phase of the trajectory. It is well known that 



Let 


•Then 



Sx (-t (38) 

and 

A(t a ) < 5 x(taJ - A 6 fcb) (-ta) ' (39) 

Therefore, since is arbitrary, 

X (ifi) = $ j (40) 


and the A S can be transformed backwards across a coast as soon as 3? 
is determined. We adopt a convention of Pitkin [ Z\ , letting X represent 
the state at {. c and y the corresponding state at 4:. Therefore, 
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The only difficult part of this is the calculation of the partials 

of 

we begin with t . 


g with respect to the state K . Taking them one at a time, 


/U t)2. To ~ AS~z- (44) 

r 0 " r 0 


and 



AS 



Also, 


<5^2, a^f - ^ aSjj <9<X 

ax TT dX aoT ax 


( 47 ) 


One can show the following: 

dS 


dS c 


n cj 

W" 5n -' ' 77 


*S, 


±Sn = 1 "v Sn + I -nS nta 
atx g t 


aS 

la 


5 = (s 3 + °< 


aS; 


aa 


0 


} (48) 

\ 


The partials of V with respect to X proceed from the fact that the coast 
time is fixed, i. e. 


9 C-fc - to) 


/* = O = 




aS, 


aS ; 


+* M 


aS. 



l ax 

ax 

ax 

l 

t S, * s 2 

a 0~ o 

Xx - 

• 


ax _ 

Therefore, 




(Yo S 0 4 - (T^j S ) 

+ ^) = -[( f " 

aS, 

+- 

a<x 

'a S 2 
r » '*x 

di 

a S3I a o( 
+■ M / — - 

aa / ax 

+ s, ar ° 

* ax 

■+S a ^ 

ax 


(49) 


(50) 


and ^ ^ is given by 

zf 


ax 


-f[Yr 0 ^Lfr e 

r Lx a<x a<x 3a/ ax ax 2 


it 

ax 


(51) 


Recalling from Eq. (48) that 


dS, 

Z&& 
a a 





1 

da Z 




' 

>' ( 52 ) 

/ 
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and defining" 


ar . 




= r n 


+ r D m aS3 


aoc 


a<X 


a <x 


(53) 


we substitute in Eq. (51) to obtain 


i± 

aX 


r 


a / d C( q 3 Q d G~o 

-t Of — y - f- Oz 


a<x ax 


ax 


ax 


(54) 


From Eq. (47), 

5 S 2 _ _ Si <?'T _ S,* dr c 

ax 


S,S 


r a^x ax 


a x 


1 °2- 
v 


d<r c aS 

t 


ax 


a<x 


a<x ax 


aP 

Hence, from Eq. (46) ™r is given by 

a-P 0 -P) ar 0 m-S, aT 


ax 


«>* ^ 


a tt- 


ax 


Tf 0 a<x ax 


f- 


+ 


rr a ax 

/j.5j Sz a$T _ a a<x 
vr 0 ax" ~ y ^ a^T ax - 


(55) 


(56) 


or 


H 

ax 


V/-r) 

l"^T" 


v 

- ps. 


ax' 


- PS: 


a n 
2 ax“ 


P f A. aS z 


a<x 


n> a cx 


aoc 
a x 


Eq. (57) can, of course, be written: 


df _ aP ar 0 aP an af acx 

3 (To ax a<x ax 


ax 


ar 0 ax 


(58) 


where 


a£_ 

an, 

5f_ 

an 


( 1 - P) 


To 


- - P'S; 




^ a t ^ A a 5* 


a# 


a<x r c a* 




(59) 



(57) 
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il- 

5 X 


Hence 


where 


i£_ 

aX 

aX 

ax 


Now we derive the partiais of with respect to X : 


- 

ax 

- M 

«?S* _ 

ax 

-mA 

a7\ 
ax " 

M 

<?£X 

acx 

ax 

(60) 

mS 2 

?r 

a fX 

+ 

m5z c 

i ar a 

- f 

^S- 

£ a^ 


r 

a oc 

ax 

V 

,j — 

r 

ax 







aS 

3 

a<x 






-/a 


ax 

(61) 

il 

. ^ 

i 

i ^ 

+ il 

a>f<j 


a<x 


- 

ax> 

a r 

e ax 

3<S 

ax + 

aox 

ax 


(62) 


IX' °'9 )Si 

Sr <’-s>& 


22L_ 

ax 


’ jT aS 3 


acx acx 

Q 

We derive now the partiais of f with respect to X • 

. ^S.C'vfjr *- r 7^) 

(n,) 1 

f ( 


s, 


rr: 


ax 


r S, 

/U — ‘ 


s 


ax 
a 

1 ;x 


a r 6 . 

r r< 


ax 


« . _ ar j 


f (r, 


aSo 
0 acx 


i 0*3- 


c ax 

, 3$ } 


( 63 ) 


o a<x 


■f M 


( 64 ) 

( 65 ) 

aS^\ a<x 


9 


a<x/ ax 


t (^o So 1 t' UY^) S/) ^ (66) 
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The only new term here is 


sS c 

d (X 



f- cX 


a Sx 

d Oi. 




Therefore, 

iL 

ax 


Q * T o 

57 


. o ar 0 zy a<X jr ji- 

- i- O, — - f y— — f- 


ax 


a* ax aT ax 


where 


a X O ,, — a *S ,• 

M ' r oS2. t- <To J— 


O ~T— < ^ 2 - 


ac< 


}(67) 


ar 

TF 

a S', 

a x 


= fc Sc i~ Cpt t r c <x) Sj 


= S 


t pS, ^ 


ax 


f- 


ac< a x 


( 68 ) 


Substituting Eqs. (67) and (68) into Eq. (65) yields 

j-f- _ mS c a ^ a S', a <x I? 5r c 

ax “ rr 0 ax' " rr 0 ToT ax ~ rj - ax 

n s>r 0 _¥ „ ar c f Ji' « 9- 3Y 

~ r « o ax r 1 ax r a<x ax' r JJ jx, (69) 


This can be rewritten 

jl£_ j_ / mS 
ax ’ y 


■h 


L_ r , o ar x 

r 2 ( ~ f p TFj 

/ / M ; ar \ „ 

Y 7 - ( TT + ^ Tr) 5 


ax x a T j(x 
ju Jk 

a c 


- ^ 

2 c>x " lhT 



£ 9 ac J_/ /x ^S, l p ar \ a<x 
r ' ax r X r 0 art a<x/ Jx\ 


( 70 ) 



c) X 
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Hence, 


where 


fashion: 


£l = 

JX 


» 


ax 


0 -F 5 To ^ -T 

a>T lx” + 77; 


a $7 0 P <5tX 

ax ^ <9o< ax 




af 

» 

p 

~ar 

5, 

/ 

ZS C " 


ar* 

f 

.17 

■ 

To 

r 


<f) 

ii 

i 

7 

[t 

t)So . 
r + 

P5a 

r 

ar 

7F“ 

P5, 


i 

i 

( AX <7> 

» 

i p 

arUr 

/A 

a^x 

7 

. 7 

' r* 

f +- 

7 7/** 

v c 


a7_ 

ao( 


-P 


> (71) 


ar 

17 


Now the partials of ^ with respect to K are derived in a similar 

(Eq. 35) 


* r- mSz 

3^ ■ 


r 


II 

ax 


r ( ax M ax*)“ ~ M ^'z) tt 


(72) 


ii = m 


ax 


r 


o ar a 5 

02 tet- r 


ax 


ax J 


(73) 


£1= 1 
ax r 




M 


a Tz 
ax 


/ 

7 


2r„ 




(74) 

a r a c< 

7x 77 


/ , * \ <3v~ a V ci a y £ atx 

7T ax Ab/ ax~ A a^x ax. 


(75) 


A- 13 






21- i_ 

a* ~ r 


(l-$) Sc 


ar c 

ax 


f C /~5 ^ S; 


<?<C 
ax + 



ar £o< 

<?0< C^X 



a<x 

ax" 


ar e 

c)X 


7 ( 0 - 3 ) 


IZ Jx~ 


(76) 


atj _ ^ a 3 ar « ^ ^3 

ax~~ ar n ax” ^ ax - *" acx ax 


(77) 


where 


13L. 

a<j 

a<^ 


) 

r 

/ 

T 


Y 


(i-j) S.-%-(li-3> Tf'^)} 

(I-3) S, - %-((l-j.>fy-/*S,) ] 


o-j) 4w- 
J a<x 


M — 

9c* 


1 f , , » \ a y <p \ a 

?'( (, 7 > JY~ m ^) jcx 


Using Eq. (42) we can define a matrix A as 


A = 


2A 

ax j 

£Y3_ 

a x j 


ax 

» 

O 

*/* 

ax. 


X, x 4 

X2. X5 
X 3 X &, 


af 

ax, 


a£_ 

ax^, 


... ^3_ 


ax. 


ax. 


+ 


f 


o 


( 


3 ° 

3 

0 3 


( 78 ) 


A- 14 





X, 

X - 4 . 


af 

3f 

af. 





or c 

a<£ 

dX 

*— 

X*. 

X 5 



ei 

2Z 


*3 



ar 0 


jc< 


ar; 


ax 



-t 

dA 

- 

JOC 


dX 


- 



■F o ’ j 

* ' 5 

O f " O 


a 


( 79 ) 


0 


Note that 

dYc 

aX~ 

' ilk 

3>X 

dot 

JX 

Defining £} to be 


X, 


*2 X, 


o o o 


r 


To r 0 r 0 

x+ x 5 x 6 X, x z x 3 
ZjuiX, ZM-Xt. ZMX 3 


we have 


8 - 


X, x 4 

X 2 Xs 
X a Xt 


To 3 

V"^ 
* Q 


3 


a y+ . . 


ay 4 


ax, 


6 s 

¥ 

* 

# 





*y t 


2X ( 


r“ 

-1 


— 


a f a P a P 





ax 


ax 


¥ 

* ¥ 


a*r* 

-t 


£2. £2 


a* 


ar e 

a^ ap< 


ax 


L 



ax 



j&x ^ JZK $■ 


J 


f o • j o 

f ■ 3 

* 

s f ' o g 


(80) 


(81) 


( 82 ) 


and therefore, the transition matrix $ defined in equation (37) is written: 



( 83 ) 
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If the transition matrix is rewritten in partitioned form as 



it can be shown that 
own elements, i. e. 


$ - 


-I 




1 

1 


1 

N 1 
— 1 
J. 

<t> 2Z 

sily found in 

-* 1 

^ T ' 

022. , 

| 

'<P,l 

! 



(84) 


(85) 


APPLICATION TO SPACE TRAJECTORIES - 
VARIATION OF PARAMETERS 


The equations developed thus far can be used to develop a very 
straight-forward Variation of Parameters integration method. The para- 
meters to be varied of course are the six initial conditions X 0 , X 0 . Again 
letting X represent the state (x 0 ,X o ) at b 0 and y represent the state at 
time t we can write 

±.-_ v il (86) 

di a-t <ay dt 

Treating the time derivatives as the sum of a two body part and a 
part due to perturbative accelerations, Eq. (86) becomes 


dx 

dt 



£*_ (i y) -hf^) 

<*y Lxdty^g vdt/p 


(87) 


Since- •>< is constant on a two body orbit the sum of the first two 
terms is zero, i. e. 
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^x\ 


aX /dy > 


( 88 } 


w 


Now, an osculating two body orbit will match position and velocity 

/dy 


ith the actual orbit, and since the first three elements of are y 
equation (87) becomes 


dx 

dt 


dX 

^7 


o 


* # ✓ # 


(89) 


wh 


ere <51 p are the perturbative accelerations. 


In terms of the notation of Eq. (85), Eq. (89) may be written. 

r -0* 

+ • * * 

K 


dx 

d t 


(90) 


Since (J;^ and 0 (l are themselves functions of X and "t and 
since <dp is also a function of x and b through the transformation given 
by Eq. (41), Eq. (90) is a well defined differential equation for X . 
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APPENDIX B 




where S here is the 6x1 vector of spherical components 



and P is the 6x1 vector of plumbline components 
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The matrix N may be partitioned into four 3x3 submatrices, 

i. e. , 


N = 


N 11 

N 12 

N 21 

N 22 


These submatrices are defined by the equations: 




3w 

s 

3w 

s 

dw 

s 

3w 

3u 

Sv' 

3u 

s 

3u 

s 

3u 
s 

3w 

3u 

3v 

Sv 

s 

CQ 

i> 

SO 

3v 

s 

3w 

du 

9v 


3w 

3w 

3w 

s 

s 

s 

3x 


3z 

3u 

9u 

du 

s 

s 

s 

3x 

3y 

3z 

3v 

3v 

3v 

8 

s 

s 

3x 

Sy 

3z 
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where 


Sw 

~aT " <^ 2 «-a 22 v-w s (d 13 cose + d 12 sin8»/(r sine) 


3w 

9 f = (n 12 v-a 32 w- Ws (d 23 oos0+d 22 sine))/(rsine) 


9w 

TT = (a 22 w-a 12 u-w s (d 33 cose + d 32 sme»/(rsine) 


dtt 
s 

dz 


(w - d U )/r 
lz s 


Su 
s 

3 x 


(u ‘ d 22 u s )/r 


du 

= <v '' d 32 u S )/r 


dv 

a# “ (d u w s ctne ' d l3 u s )/r 


3v 

s 

ay 

av 
s 

dz 


(d 21 w s Ctne - d 23 u s )/r 


< d 31 W s CtnS ' d 33 U s )/r 
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APPENDIX C 



FIRST PARTIAL DERIVATIVES OF 
GRAVITATIONAL ACCELERATION WITH 
RESPECT TO PLUMBLINE POSITION COORDINATES 


The matrix J is defined to be the matrix of first partial derivatives 
of the gravitational acceleration vector in the phimbline system with 
respect to the plumbline position coordinates. This matrix is used in the 
gravity related terms of the adjoint (Euler-Lagrange) equations. 


\ 

3 y 

9g x 

dz 


g xx 

g xy 

g 

& xz 

3 g y 

3 g y 

Sg y 





9x 

3 y 

3x 


g yx 

g yy 

g 

yz 

if*. 

Sx 

3g z 

3 y 

3z 


g zx 

g zy 

g 

& zz 


J = G U I + 


X 


12 


22 


32 


G 22 

G 23 


X 

tSJ 

1 

°3 2 

G 33 

• 

_ a i2 

a a 

22 32_ 


is defined in Section 3.1, I is a 3x3 identity matrix, the 

a are elements of the A matrix, and 
1 ] 


3 © 
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i ^ G .. , Q 9 G 1 1 

p JL , 1J. ctn0 11. 

<?» ' *\ ._ „ "N n • 


22 r ' 9r 


90 


At /R \ 2 
f p E 

2 L 11 3 

r r 


/^ e \^ 9 20 2 /® e Y^ 2 

(Cj{ (| - ~ cos 1) + H ( r ) (4" 14 cos 0) cos 0 


(y- <12-24cos^0) cos^0))] 


+ DJ 


G„„ ~ G 


23 3 2 r sin.0 90 


1 5G 11 _ 1 , 5G TO ctn0 5G TO . 


r 9r 


90 




( R \ 2 ,• n \ - „ 

~ J cos0 - H (3 - 21 cos 0) 


,R x 3 


/R \4 


+ DJ (~~j (12 "36 cos 2 0) cos0] 


G. 


1 sg to 


33 r sin0 9 0 




|[2CJ^) +6h(-^) cos0 +-y Dj(”) (1 -7 cos 2 0)] 


The fact that J is symmetric can be anticipated, since J is also 
the matrix of second partial derivatives of the gravitational potential 
function U(r, 0) with respect to the plumbline position coordinates. 


02 



to 


with 


since 



In the event that a spherical earth is being simulated J reduces 


G u 1 + 


G 22 [ x y *] 



G 23 ' G 32 ' G 33 = 
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APPENDIX D 


RELATION OF C A AND C N TO C L AND C D 


The aerodynamic formulation in Section 5 is predicated upon the 
assumption that the lift and drag coefficients may be approximated by 


and that 


C L * C L “ 
a 


c d' c d +c d “ 

o a 


a & sina 


a fid 2(1 - cosa) 


Using these assumptions and become 


C = C sina 

Xj J-j 

a 


C D = C D + 2C d (1 - cosa) 
o a 


Since 


cosa - sina 


= C T cosa + CL sina 
N L D 
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there results 


C = (C + 2 C ) cos a + (2 C - C ) sin 2 # - 2 C 
ADD JJ J-j J-) 

o a a a a 


C = (C + 2 C ) sin# - (2C - C ) sin# cos# 

N D D D L 

o a a a 


Therefore the aerodynamic coefficients a, b, c of Section 5 are 


a = C D + 2 C d 
o a 


b ■ 2C D - C L 
a a 


c = -2C 


D 


D. 1 REDUCTION OF DATA 

Given C (a, M) and C (a, M) it is easy to construct 

A 

C^a, M) = cos# - sin# 

C~(a, M) = C. cos# + Q, sin# 

D A N 

therefore, C and C may be thought of as perfectly general. 
J-J JJ 
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It is not unreasonable to expect to be able to approximate a given 
set of C (a, M), C (a, M) data as 

J-j JD 


C T = C T +C T a* 

j-i j-i S-i 

o ' a. 


C D - °D + C D, C L, 
o L 


dC 


where 


D 


D 


L dC, 


and where a* denotes the reference or data-base angle of attack. 


Since the angle of attack is the angle between some reference 
direction and the velocity vector, define a new reference direction by 
the angle of attack, a, related to the one on which the data is based by 


a = a* + 6 


where 6 is a bias angle. 


In terms of a , C is 

J-i 


C L ■ C L 


In order that C = 0 when a = 0 

JLj 


+ C L (a -6) 
o a 

it is necessary that 


6 
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In general 6 will not be exactly constant with respect to Mach 
number, however, a good least squares constant bias may be calculated 
by choosing a set of Mach numbers, M.., and then forming 

6 =. EC t (M.) • C T (M.)/SC t (M.) 

. 1 JLi 1 . _Lr 1 


In terms of a, becomes 


C D * C D +C D < C L +C L 
o L o a 


- C D +C D «°L - C L 6)+C L a) 
o L o a. ol 


Recall that 6 was chosen such that 


C L - C L 6 - ° 
o OL 


Therefore to good approximation 


^ ^ 2 2 

C D " °D C D t °L “ 
o La 


or 


S = C D +C D “ 

o a 


where 


C D 3 C D t C L 
a La 


which is consistent with the original assumptions about the form of C 


D* 
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APPENDIX E 



RADIUS OF APOGEE AND RADIUS OF PERIGEE CONSTRAINTS 


The radius of apogee and radius of perigee constraints are treated 

in the APPG package which includes subroutines APPG, ADDR, AGE($, 

F0RKM and SYSDER. The approach used is to compute the radii of 

apogee and perigee r and r , respectively, and the associated adjoint 

a p 

vectors X . and X, through a numerical integration of the orbital 

r* ' 1C 

x a p 

equations of motion and the associated adjoint equations. The orbital 
Euler-Lagrange equations are 


x = !(x) 



5H 

5x 


T 

) 


where H is the variational Hamiltonian. These equations take the form 


w 

= 

G 

X 1 = 




X 

1 

4 

u 

— 

G 

X n = 

"X c 



y 

2 

5 

♦ 

V 

zr 

G 

X 0 = 




z 

3 

6 

X 

zz 

w 

X, = 

-X.G - X G ~X 0 G 




4 

1 xx 2 yx 3 zx 

♦ 

y 

- 

u 

X, = 

“X. G - X 0 G "X^G 




5 

1 xy 2 yy 3 zy 

z 

- 

V 

X c = 

-X. G -X 0 G - X 0 G 




6 

1 xz 2 yz 3 zz 
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where (G , G , G ) is the gradient and G , G , etc. are the second 
x y z s xx yy 

partials of the gravity potential. Boundary conditions are x(t Q ) the state 

at injection given and 



dr 

dx t = t 
— a 



dr 

dx t = t 

P 


where t and t represent the times at which the orbit passes through 
a p 

its apogee and perigee, respectively.- 


The procedure used in the numerical integration of the orbital 
equations is to first integrate the system equations forward in time from 
x(t Q ) to some time t^ at which time r passes through its first extremum. 
The stopping condition used is 




r • 


r 


t=t 


[xw+yu + zv]^_^. 

1 


1 


0 


The adjoint equations are then integrated backward in time from t^ to 

t from the terminal conditions 
o 


T 

x. <v 


dr 

dx 


t = t„ 
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Define X, as the solution at time t of this backward integration. 
r l 

The second extremum is found by integrating the system equations 

backward in time from x{t ) to some time t at which time r passes through 

o z 

its first extremum. The adjoint equations are then integrated forward in 
time from the initial conditions 




Sr 

Sx 


t = t 


2 


Let X, be the solution of this forward integration of the adjoint equations 
Y r 

2 

at time t . Then a simple comparison of r(t 1 ) and r(t ) can be used to 

v i u 

set r , r , X, and X, .If r(t ) is greater than r(t ) 
a P ~V r ~r r z 1 

a p 


r = r(t ) X , = X. 

p i ~i> -i> r 

p i 


r a ' r(t 2 } b ~ b 

r a r 2 


Otherwise 


r P ' r( V b ’ b 

r p r 2 


r a = r(t l ) b = b 

r a r l 
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The advantage of using this approach to determine r^ r^, 


and X, is that it does not require analytic expressions for r^, r , etc. , 
^r 


P 

in terms of the state at injection. In general derivation of analytic express- 
ions for r , r , etc. , for an oblate gravity field is quite difficult and some 
a P 

sort of approximation is generally used. By computing these quantities 
through a numerical integration of the orbit equations, this approximation 
is avoided. 


The only limit on the accuracy for this approach is the numerical 
integration error. Numerical integration of the orbit equations in APPG 
employs a fourth order Runge Kutta differential equation solver with a 
fixed step size of 20 sec for the system equations and 40 sec for the 
adjoint equations. Numerical accuracy appears to be adequate but can 
be improved if desired by decreasing the step size DT. 
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APPENDIX F 

THE IMPACT POINT PENALTY 


The impact trace of a launch vehicle will be forced to avoid certain 
regions by utilizing an impact trace penalty function. A new state variable 
is therefore defined to be 


' n+1 


k 

= E 
i-1 


C. e 

l 


T 

-Z. B.Z. 

ill 


where CL are weighting factors. Eh is a 2x2 symmetric matrix associated 
with each point ( G^ 9 ) that shapes the exponential function, and 


Z. 

l 



By properly choosing C., 0.^ and 0^ the impact trace penalty function 

can be designed to cover forbidden zones in the 0, 0, (latitude, longitude) 

space. By constrainting the termal value of x the impact trace can be 

n+1 

made to avoid the "peaks" in x . An efficient algorithm for calculating 
the latitude and longitude at impact which is consistent with the partial 
derivatives needed for the adjoints is developed below. 


In order to avoid iterating on the change in eccentric anomaly 
necessary to give the proper coast time, it will be computed analytically 
and then used in the Goodyear coast equations. 

Figure 9 -shows a schematic of the coast to impact. 
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FIGURE 9 SCHEMATIC OF COAST TO IMPACT 
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Obviously, 



and 


„ -1 ,a - r x 

E = cos ( —~i 
n ae 


a- R 

e t ° 2 * ' cos < ~7T ) 


It is assumed here that only coasts initiating with positive flight path 
angles are of interest. 

The charge in eccentric anomaly is 


AE = 2ir ~ (© 1 + 6 ) 


where 


a . a - R 

A -1 , e. 

= cos ( ) 

1 ae 


„ A -1 ,a ~ r x 
6 2 ' C0S 


Thus 


cosAE = cos(0 1 + 0 2 ) 
sinAE = -sin(6^ + e^) 

Since both 9 and 6 will be calculated to be less than pi, the sine will 

X & 

be positive for both angles; therefore we can calculate 
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sinAE = -[cos 9^ 


v 2 yi 


with. 


cos 9 . 


a- R 
e 

ae 


cos 0 2 


a - r 
ae 


The correspondence between AE and Goodyear's universal variable, 
4 !), is known to be 

ip = (a/Ai) 1/2 AE 

Also, 

S = C = cosAE 
o o 

S t = *C = (a/ju)^ 2 sinAE 

5 2 =0 2 C 2 = (a/juAl“S o ) 

53 =^ 3 C 3 = (a//i)(^-S 1 ) 

,2 

s 4 =0 c 4 = (a/^Kfr - s 2 ) 

3 

S 5 ^ 5 C 5 = (a/jaXfi- - S 3 ) 
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Then defining 



a - (r- v). ... „ 
o initial 

along with 

F = 1 - juS 2 /r 
g = rS 1+ a o S 2 

the position coordinates at impact are 


h\ 


I X ^ 


M 

* 

= f 

y 

+ g 

u 

1 z 1 

1 

z I 


\ v l 

\ p/ 


\ / 



The time to impact is 

At = rS +a S„ -1- juS 
1 o 2 3 

Using the A matrix of Section 2, 


l x) 

i • 

T 

fx \ 

J 

= A 

y 

\ z ! 

1 

\\i 
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the latitude at impact is 


e j> = sin X (y/R > 

P 

and the earth relative longitude is 

0 = tan _1 (x/z) - O (At + DTZ + T) - AL0NG0 

Xt 6 

P 

Using these values, i? n+ ^ can be calculated. 

During the adjoint run, the partial derivatives of with 

respect to the state along the trajectory must be calculated. 


The calculation proceeds as follows: 


dx 


T 

sr = ' EC i e iZi(2b iii ( %- e i> + b ii2 ( V 0 i ) 


. . T_ 

ox . “z. B.z. 

TT S ' EC i C 1 11<b i 1 2<V 9 i > + 2b i22 ( V (S i ,) 

P 


dr = 0 = 


5x + r dt 
Sx^ p f 


_-l dr j, 
dt„ = v - — ox 
f r bx p 
P 
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where 


r = r 


r /R 
P e 



" 3x ' 3x 

E • £ 

_ 3x 3v_ 



30 30 

d0 = r— ^ fix + r— ^ x dt 
p 3x p 3x p x 
P P 


d0_ l 

d0 = - — fix + — — ^ x dt „ 

p 3x P 1 3x p e f 
P \ P 1 


which can be written 


30 r 


d0 = 


p L 


i - 


1 • dr I e 

T x r l fix 

r p 3x j p 

p j i' 


d0 = 
P 


J^£ 

I - ~ 

3r 

x 

n 

+ - 4 - 

3r 1 

) 3x 

r 

p 3x 

r 

3x f 

l P 


P 


P J 


fix 


• B r 

Note that x - — is a 3x3 outer product 
p 3x 



R cos 0 
e p 


A“ 


~x 

2 2 
x + z 
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or 


dx 


n+1 


'3x 


n+1 


30 


n+1 , dX n+l ,, 

dx = — — — d0 + — — d0 

n+i p 90 p 


.p x +z 


3x 


P 


n+1 


P 


_ 3x 


n+1 


2, 2 30 R cos 6 


P 


P 


30 


x 


2 , 2 
x + z 


T 

a r 


where 


, d^rr+l 3r . 

30^ r 3x_ f“ 6x p 


■p A T 1 • 

1 - I„ ~ — x - — 
3 r p 3x 


Some simplifications in the numerical work required can be made. 
Defining 


A SX .1 

A n+1 
a = 


30 2,2 

p x +z 


A 11 

, A n+1 
b = — — - * 


3x R cos 0 

pep 


This gives 


3x 


n+1 


9x 


P 


3x n+l 

= <PR : pr 2 p R 3 )r + 


Q e pT 
r *v p 


P P P 
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where 


Then defining 


gives 


PH 1 = < Za u' Xa 13 )a + ba 12 

PR 2 = (Za 2! ' Xa 23 )a + ba 22 

PB 3 = (Za 31 -x a33 )a + ba s2 


c 4 (PR i + PR y + PR i )/(r • v ) 
1 P 2 * p 3 p PP 


a 3x , 

. A n+1 
d = 


a 


30 , 


r . v 
P P 


3x 


n+1 


3x 


P 


(PR, +(d~c)x PR 0 + (d-c)y PR + (d-c)z ) 
1 p 2 •'p 3 p 


Finally 



3x , . 
n+1 

3x 

n+1 

3x 

P 

t 

t dx 
t E 

Sx 

5v 

3x 

P 

3x 

| 3v 


where the partials of state at impact with respect to state along trajectory- 
will be calculated using Goodyear’s partials. 
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The velocity components at impact are calculated from 
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APPENDIX G 


INPUT DESCRIPTION AND EXAMPLE PROBLEM 


The user of the LIFTING ROBOT program will find it helpful to 
sketch a thrust profile before setting up input for the problem he wishes 
simulated. Sketched in Fig. 9 is an 8 thrust event representation of a three 
stage Saturn V thrust history. Vertical lines and horizontal lines will 
be referred to as "pickets" and "spaces", respectively. The "picket" 
numbers in Fig. 10 are circled. Note, there is always one more picket 
than spaces. A thrust event must be defined every time there is a dis- 
continuity in the total thrust. Dashed vertical lines represent miscella- 
neous weight drops. Spaces are thrust duration times and are denoted 
TAUT. The elapsed time between the Jth miscellaneous weight drop 
event (dashed vertical lines) and some thrust event picket is denoted 
TAUW(J). The particular thrust event picket to use is denoted N0WD(J). 
LIFTING ROBOT drops the atmosphere at the IWDCHI th miscellaneous 
weight drop event. Therefore, a miscellaneous weight must b e dropped 
where the atmosphere is to end even if it is a zero (0) weight drop. 

Min-H optimization begins at either t or t~ whichever occurs first. 

X ^ 


The LIFTING ROBOT program controls vehicle flight by looking 

up x and X as a function of time out of control tables. The Min-H 
3? if 

steepest ascent process adjust these tabular points until they take on 

optimal values. A "control table" consists therefore of three tabular 

arrays: time, x * X • LIFTING ROBOT contains four control tables, 

3? y 

each containing a maximum of 49 points. In order to provide generality 
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FIG. 10. SATURN E THRUST PROFILE 
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for the user, the Jth control table begins at the NBGCT(J) th picket, ends 
at the NENDCT(J) th picket and has a maximum of NP(J) < 49 points. NP 
should be odd for all tables in use and zero for all others. Control tables 
should not extend over coasts or over an intermediate constraint point. 

If Min-H is to begin in the middle of a thrust event, NBGCT{1) should be 
set to the picket at the beginning of the thrust event. 
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NAMELIST INPUT DESCRIPTION 


INPUT ,, 
SYMBOL “ 

SIZE 

EXPLANATION 

PRESET 

VALUE 

UNITS 

HEAD 

(15) 

Identification for print out (60 
characters) 



TZERQ) 


Initial time 


sec 

TLIFT 


End Of lift-off; beginning of tilt 


sec 

TTILT 


End of tilt 


sec 

TCHFRZ 


Begin Min-H 


sec 

DTZ 


Time from GRR to Lift-off 


sec 

FLBS 

(15) 

Thrust per engine/ thrust event 

0 . 

lbs 

TNE 

(4, 15) 

Number of engines /thrust event 

0 . 



Four numbers for each thrust 
event: the number of inboard 
engines, their cant angle (deg), 
the number of outboard engines, 
their cant angle, (deg). 


WD0T (15) 

Flow rate per engine /thrust 
event 

0 . 

lbs /sec 

CW0T (15) 

Critical flow rate per engine/ 
thrust event 


lbs/sec 

'$ INPUT in Col. 2. All Data begins in Col. 2 
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INPUT 

SYMBOL 

SIZE 


PRESET 

VALUE 

UNITS 

WDLBS 

(15) 

Weight dropped during a weight 
drop event 

0 . 

lbs 

WTJET 

(15) 

Jettison weight/thrust event 

0 . 

lbs 

AE 

i 

(15) 

Engine exit area/ thrust 
event 

0 . 

2 

m 

S 

(15) 

Aerodynamic reference 
area /thrust event 

0 . 

2 

m 

TAUT 

(15) 

Thrust event duration time/ 
thrust event 

0 . 

sec 

TAUW 

(15) 

Elapsed time between a 
thrust event and a weight 
drop event 

0 . 

sec 

N0WD 

(15) 

Denotes a picket number from 
which TAUW is defined 

0 


N0EVNT 

(5) 

The total number of thrust 
events which comprise a 
stage 

0 


PRINT 

(15) 

Print time increment /thrust 
event 

10. 

sec 

STEP 

(15) 

Integration step- size increment 
for forward run/thrust event 

8. 

sec 
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INPUT 

SYMBOL 

SIZE 


PRESET 

VALUE 

UNITS 

BSTEP 

(15) 

Integration step- size 
increment for backward 
run /thrust event 

16. 

sec 

AZ 


Launch azimuth 

90. 

deg 

LAT 


Initial geodetic latitude 

28. 531855 

deg 

XJEXT 


= 1. if maximizing payoff 
"-1. if minimizing payoff 

1 . 


CASE 


Case number 

1 . 


DP2 


Decimal fraction of constraint 
error to remove in first 
iteration 

. 5 


QY 


Decimal fraction of H to 

a 

remover per iteration 

. 8 


CHIDT 


« 

X for tilt -over during first 
stage pitch 

. 1 

deg/sec 

W01 


Lift-off weight at TZER0 

0 . 

lbs 

DELVG 


AV for geometry reserves 

0 . 

m/sec. 

DELVP 


AV for performance 

0 . 

m/ sec. 


reserves 
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INPUT 

SYMBOL 


SIZE 


EXPLANATION 


PRESET 

VALUE 


UNITS 


WPM 

Maximum critical propellant 
in stage from which 
performance reserves are 
taken 

0 . 

lbs. 

TCHIR 

Time of chi roll initiation 
(for report tables) 


sec 

CHRD0T 

Roll rate (for report tables) 


deg/ sec 

FAZ 

Azimuth at which Fin 1 
points (for report tables) 


deg 

AL0NG0 

Longitude of the launch site 
(measured positive west) 

80. 5649528 

deg 

EU 

Upper error bound in forward 
integration 

1. E-5 


BEU 

Upper error bound in 
backward integration 

2, E“5 


AYL 

Used for error check in 
forward integration 

2 ,. E-3 


BYL 

Used for error check in 
backward integration 

4. E-5 


HMN 

Minimum step-size for 
forward integration 

. 25 

sec 
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INPUT 

SYMBOL 

SIZE 

EXPLANATION 

PRESET 

VALUE 

UNITS 

BHMN 


Minimum step-size for 
backward integration 

. 50 

sec 

CMUE 


Gravitational constant 

3. 986032E14 

3, : 

m /sec 

0MEGA 


Angular rotational 
velocity of earth 

' 7.2921158E-5 

rad/sec 

CJ 


First coefficient in 
gravitational expansion 

1. 62345E-3 


H 


Second coefficient in 
gravitational expansion 

-5.75E-06 


DJ 


Third coefficient in 
gravitational expansion 

7. 875E-06 


FLAT 


Flattening of Fischer 
ellipsoid 

1/298. 3 


RE 


Equatorial radius 

6378165. 0 

m 

GZER0 


Relates mass to weight 

9. 80665 

m/sec 

J0RB 


= 1 if spherical earth 

= 0 if oblate earth 

1 

0 


JUMP 


Jump start at this picket 
number if JUMP > 1 

1 
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INPUT 

SYMBOL 

XWDCHI 

KIND 

KWTA 

NMAX 
NT ABLE 

NVRST 
I PR 


SIZE EXPLANATION 


PRESET 

VALUE UNITS 


The number of the 
weight drop event where 
atmosphere is dropped 

Type of integration used: 

= 1 for variable step 
size Adams- Moulton 
= 2 for Runge Kutta 
=3 for fixed step Adams 

=2 if only optimized 
=3 if X and \ optimized 

p y 

Total number of iterations 

= 1 if output tables are 
wanted for publication 

Intermediate constraints 
imposed at termination of 
this thrust event. Must be 
zero if no intermediate 
constraints wanted. 

Thrust event from which 
performance reserves are 
taken. <IPR must be zero 
if no performance reserves 
are wanted). If IPR^O, WPM 
and CWD0T must be input. 


1 


3 


2 

0 

0 

0 


0 
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INPUT 

SYMBOL 

LAST 

KRDER 

KINDB 

KDERB 

IMP 

JTHR 

KCDPHI 


SIZE EXPLANATION 


PRESET 

VALUE UNITS 


=0 if only one case is run; 0 

= 1 if more cases are run 

Order of differences in 3 

integration package for 
forward run. 

Type of integration used in 3 

backward' run (See 

JRBETC+3) 

Order of differences in 3 

integration package used 
for backward run. 

Jettison weight of this thrust 0 

event will be integrated to 
impact. (Cannot be the last 
thrust event). 


(15) =0 if input thrust and flowrate 0 

used 

= 1 if thrust and flowrate 
found in ATTRAC 
= -l if thrust and delta 

weight found in ATTRAC 

(10) Terminal function codes. 

(Code in KCDPHI(l) is 
payoff) 
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INPUT 

SYMBOL 

SIZE 


PRESET 

VALUE 

UNITS 

PSIREQ 

(10) 

Constraint values desired 
at terminal point. (Value 
in PSIREQ(l) is constraint 
for code in KCDPHI(2), etc. ) 



KCDRES 

(6) 

Intermediate constraint 
function codes 



PSIRST 

(6) 

Constraint values desired 
at restart point. 



RDB 


Control parameter switches 



KDB{1) 

KDB{2) 

9 


Insert 1 to optimize TAUT(l) 
Insert 1 to optimize TAUTX2) 

4 

0 

0 


• 

• 

KDB(16) 

KDB{17) 

KDB{18) 


m 

Insert 1 to optimize W01 
Insert 1 to optimize CHIDT 
Insert 1 to optimize AZ 



KDT 

(17) 

Companion vector to KDB. 
Contains in corresponding 
locations the number of the 
thrust event from the present 
one which is to be altered in 
order to hold tank limit. 

0 
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INPUT 

SYMBOL 

SIZE 

EXPLANATION 

PRESET 

VALUE 

UNITS 

WIBT 

(17) 

Used to speed up or slow 
down convergence of one 
parameter relative to 
another, 1st element of 
WIBT goes with 1st active 
parameter, 2nd with 2nd 
active parameter, etc. 

1. 


NBGCT 

(4) 

Jth contrdl table begins at 
NBGCT(J)th picket 



NENDCT 

(4) 

Jth control table ends at 
NENDCT(J)th picket 



NP 

(4) 

The number of points in a 
control table (Must be an 
odd number of points, ) 

0 pts 




CONTROL TABLES 



TTBL(l) 

TTBL(5I) 

TTBL(lOl). 

TTBL(151) 

49 pts. 

If 

IT 

TT 

1st time table (real time from 
TZER0) 

2nd time table (real time from 
TZER0) 

3rd time table (real time from 
TZER0) 

4th time table (real time from 
TZER0) 


sec 
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INPUT 



SYMBOL 

SIZE 

EXPLANATION 


PRESET 

VALUE UNITS 


CPTBL(l) 
CPTBL(5) 
CPTBL(lOl) 
CPTBL(1 51) 

CYTBL(l) 

CYTBL(51) 

CYTBL(lOl) 

CYTBL(151) 

VIV 


RTASC 

DECL 

VELTGT 

GAMTGT 


(49 pts. ) 1st v table 

P 

" 2nd X table 

P 

" 3rd\ table 

P 

" 4th x table 

*P 

(49 pts. ) lstx table 

11 2nd * table 

y • 

3rd X table 

y 

" 4th x table 

y 

(8) Vector of initial conditions 

for a jump start. 

If VIV(7)=0, input: x, y, z, x, y, z(z, x, y, z, x, y Apollo 13) 

If VIV(7)=2, input: V T , y, r, A , Lat., Node 

X z 


Right ascension of 

outgoing asymptote 0 

Declination of 0 

outgoing asymptote 

Rendezvous Target 0 

.velocity at node 

Rendezvous Target 0 


Flight Path angle at 
Node 


rad 


rad 


deg 

deg 

m/sec 

deg 
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INPUT 

SYMBOL 

SIZE 

EXPLANATION 

PRESET 

VALUE 

UNITS 

RTGT 


Rendezvous Target 
Radius at node 

0 

m 

INCTGT 


Rendezvous Target 
Inclination 

0 

deg 

LATTGT 


Rendezvous Target 
Launch site Latitude 
Displacement 
(Ignored) 

0 

deg 

BTATGT 


Rendezvous Target 
Position Phase Angle 

0 

deg 

NCOST1 


Thrust Event for 1st 
Analytic Coast 

0 


NCOST2 


Thrust Event for 2nd 
Analytic Coast 

0 


IAA 


= 1 calculates launch azimuth 
for rendizvous 
t* 1 uses input A 

z 

0 


IPHIT 


= 1 calculates LATTGT 
/I uses input LATTGT 

0 


L PRINT 


-0 Ignored 

= 1 Suppress x table print-out 
=2 Suppress trajectory 
print out 

0 
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INPUT 

SYMBOL 

SIZE 

EXPLANATION 

PRESET 

VALUE 

UNITS 

IC0NT 


No. of times complete 
trajectory integration 
is carried out after 
convergence 

1 


IAEOK 


Used in ATTRAC for 
decaying exit areas 

0 


L0NGPR 


Number of complete 
trajectory printouts 
after convergence 
Cannot be greater than 
IC0NT 

1 


TH1 

(10) 

Impact point ellipse 
centers-Tatitude 

0 

deg 

PHI 

(10) 

Impact point ellipse 
centers -longitude 

0 

deg 

R0TA 

(10) 

Impact point ellipse 
rotation angle (pre- 
counterclockwise) 

0 

deg 

LATWTH 

(10) 

Impact point ellipse 

latitude width 

i 

0 

deg 

L0NWTH 

(10) 

Impact point ellipse 
longitude width 

0 

deg 
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INPUT 

SYMBOL 

SIZE 

EXPLANATION 

PRESET 

VALUE 

UNITS 

N0I 


No. of impact point 
ellipses 

1 


IPCNST 


= 0 if no impact option 
= 1 if impact penalty 

integration is desired 

0 


MSWCH 

(15) 

= 1 if staging on time 
= "1 if staging on fuel 

0 


GLIMG 

(15) 

Axial "g" limit /thrust 
event 

0 

g's 

FUELG 

(15) 

Fuel consumed in thrust 
event where MSWCH =-l 

0 

lbs 
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$ INPUT 2 


TITLE - 
0FFICE - 
DATE - 
NCASE - 
SRID 
$END 


48 Columns of BCD information 
12 columns of BCD information 
12 columns of BCD information 
Fixed point case number; should be < 1000 
360 columns of BCD information 


The alignment of the codes and constraints is 


KCDPHI 

PSIREQ 

KCDRES 

PSIRST 


payoff code 

1st constraint code 

2nd constraint code 


1st constraint value 

. 

2nd constraint value 


1st constraint code 

2nd constraint code 

1st constraint value 

2nd constraint value 


etc. 




EXAMPLE PROBLEM 



Maximize payload into perigee of a 50-100 nm orbit having an 
inclination of 55 degrees. Launch at a 38 degree azimuth from Cape 
Kennedy over an oblate earth using a two- stage space shuttle vehicle. 
Both stages must be throttled to maintain a 3 "g" axial acceleration limit 
and are staged on fuel depletion. Controls to be optimized are: lift-off 
weight and tilt-over \ as well as Xp a nd X y a fter 140 sec. 
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Data for this problem is given below: 


Thrust Event 

1 

2 

3 

Thrust /engine (lb) 

520000. 

0 . 

597000. 

Flow rate /engine 

1298. 4 

0 . 

1300. 6E 

(lb/ sec) 

Jettison weight (lb) 

700000. 

0 . 

0 . 

Number of engines 

12 

0 

2 

/ 2. 

Engine exit area (m ) 

2. 1869 

0 . 

0 . 

Aerodynamic ref. area 

929. 

929. 

0 . 

(m ) 

Burn times*(sec) 

205. 

10. 

297. 

Integration step (sec) 

1 . 

2. 

4. 

Print Interval (sec) 

20. 

20. 

20. 


Lift-off time = 0 

Begins tilt at 8 secs 

End \ tilt at 30 secs 

Begin optimization at 140 secs. 

Group thrust events as follows: 

1st stage - 2 thrust events; 2nd stage - 1 thrust event 
Start 1st control table at picket 1 and end at picket 2 (25 points) 
Start 2nd control table at picket 2 and end at picket 3 (5 points) 
Start 3rd control table at picket 3 and end at picket 4 (35 points) 



• Just guesses since staging is on fuel 
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3rd table 


|d^ 


Estimate Y and v in the three tables to be: 
P 7 


1st table 


2nd table 


Y from 1. 14 to 1. 22 from 1. 34 to 1. 37 from 

v from 0 to-. 0015 from 002 to 0025 from 

A y 


1. 23 to 1.74 
0017 to 0018 


Estimate starting tilt over \ to be . 155 deg/sec 
Estimate launch weight to be 47 88230 lbs 
Injection conditions are: 

Vel. = 7876. 4195 m/sec 
Gam = 0 

R = 6470762. m 
Incl. - 55° 


There is a 10 sec coast between the first and second thrust 

events. 


Report tables are to be output. 

This problem converges in 6 iterations. All constraints are met 
within a small tolerance and the max payload is 352840. 8 lbs. 

A listing of the input cards and aero data for this problem are on the 
following page. 
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lifting-robot namelist input data 


|DRCj 

% INPUT 

TZEKO=C.» TLIFT=8.. TTILT=30.* TCHFRZ=140 . * DTZ=0.* 

FL13S=520000.»0. *597000. * 

T NE=12 . *Q.*O.»O.*O.*0.*O.*O.*2.* 

WUOT = 1298.H»0. * 1300 .6536 * 

WTJeT=7U00u0. * 

AE=:2.1869* S=929.*929.* 

TAUT=205. *10. *297. * 

T AUfl=12 • * N0WU=2 * 

NOE VNT=2 * 1 * 

PRlNT=20. *20. *20. * 

SYEP=1 . *2. *4. * BS1EP=2. *4. *8. * 

AZ=58.* LAT=28.531* 

XJEXT=1.» CASE=1.» DP2=.5* QY=,75* 

CHIuT=. 155* 
l v01 = .478823E7* 

AL0uG0=80.o3S4* 

K lNu-1 * KWTA=3* NMAX — 1 0 * 

NIAULE=X» 

KCLkH 1=1*2*3*4.10* 

PblKEQ=7d7o. 4195*0. *6470762. .S5. . 

KD6 (16 ) =1 * 1 * 

.v 16 1=2 • * « i * 

Nt5GCT=i.2*3. NENDC I =2 *3*4* 

NP=25 * 5 * 35 * 

TT6l=140, *200. . CP l'BL=l« 14.1.22* L YTBL=0 . . - . QU 15 . 

TTBL(51)=20U. *210. . CPTBL ( 51 ) =1 . 34 . 1 . 37 . CYT0L ( 51) =- . U02 . -. 0025 * 

TTBL ( 101 ) -21 0 . .4 75. . CP TBL ( 101 ) =1 . 23 * 1 . 74 * CYTBL ( 101 ) =-. 0017*-. 001 8 * 

ICOi.T=l* LONGPR-1 * 

M S ft C H = — 1 . 1 . - 1 * 

GLlFiG— 3. . 3 . * 3 « . 

FUELG=3U53849.6*0. *681542.5. 

SENi^ 

SINPUT2 
NCASE=1 * 

SEND 


SUBROUTINE CACN(M.L) 

AERO DATA — A(M)*B<M) AND C(m> — CONF 19 BOOSTER HCP ORBITlr 

REAL M.MT 

COiMMOn/AKODTA/ACF.BCF.CCF * DAuM * D8UM * JjCDM 
CuMMOim/TABLK/K » MMM (11) 

DIMENSION MT (20) * ATB (20) *RTb(20) *CTu(2G) 

DATA N/ 12/ 

DATA ( MTU) *1=1, 12) /O. . .5* .85.1. .1.5.2. ,3. *4. *5. *6. .7. *20./ 

DATA ( ATB ( I ) * 1 = 1 . 12 ) /2 . 89t» * 4 . ci52 * 8 . 422 » 7 .917 . 4 . 781 * 4 . 853 . 

1 4.38o*3.33U*2.80y .2.651 >2.643*2.643/ 

DATA ( bTB ( I ) * 1 = 1* 12)/. 2754, l.o92*4.598 *4. 04 3*1.948*2.563* 

1 2. 7ol* 2. 069* 1.734*1.644*1.644 *1.644/ 

DATA (CTBl I > * 1 = 1 .12) /-2. 796* -4.780 *-8.322 *-7.767 *-4.641* -4. 740* 

1 -4.308 *-3.2720 *-2.765 *-2.619 *-2.619 *-2. 619/ 

CALL SPLINE (L*K*N* K * MT . ACF . ATB . BCF » 6TB * CCF * CTR * DAOM* OBDM . DCDn ) 

RE TURN 

END 
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